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FOREWORD 


The  increase  in  oil  and  gas  production  has  led  inescapably 
to  the  necessity  of  increasing  the  accuracy  of  calculations  made 
in  designing  methods  for  the  development  of  oil  and  gas  deposits 
and  of  increasing  the  reliability  of  prospecting  hole  tests.  This 
in  turn  requires  a deeper  and  more  detailed  description  of  the 
processes  and  phenomena  which  take  place  during  the  operation  of 
oil  and  gas  wells.  Until  recently,  complications  connected  with 
this  served  as  an  obstacle  to  a more  precise  and  rigorous  formulation 
of  the  appropriate  problems.  With  the  advent  of  computers  and 
the  development  of  numerical  methods  of  integration  of  the  equations 
of  mathematical  physics,  this  obstacle  may  to  a significant  degree 
be  overcome. 

In  the  present  work  the  authors  primarily  strive  to  illuminate 
several  new  questions  involving  the  nonisothermal  movement  of  fluid 
and  gas  in  pipes,  namely,  the  processes  of  the  thermal  interaction 
of  oil  and  gas  wells  with  frozen  rock.  In  addition,  the  authors  wish 
to  generalize  and  systematize  several  applied  aspects  of  pipe 
hydraulics. 

The  work  was  carried  out  at  the  I.M.  Gubkin  Moscow  Institute 
of  Petrochemical  and  Gas  Industry  and  at  the  Laboratory  of  the  Mechanics 
Continua  and  Dispersed  Media  of  the  Institute  of  Physical  and  Technical 
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Problems  of  the  North,  affiliated  with  the  Yakutsk  Branch  of  the 
SO  AN  SSSR. 

The  authors  express  their  thanks  to  their  senior  scientific 
coworker  at  the  Department  of  Petroleum  & Gas  and  Subterranean 
Hydromechanics,  V.I.  Maron,  for  his  great  assistance  in  preparing  the 
material  of  this  book,  and  also  to  the  coworkers  of  the  Laboratory  of 
Continuum  and  Dispersed  Mechanics  who  assisted  and  took  part  in  the 
investigations,  especially  L.K.  Stepanets,  who  devoted  much  effort  in  the 
layout  of  this  book  and  in  preparing  the  manuscript  for  printing. 

E.A.  Bondarev  wrote  Section  1 of  Chapter  1 and  Sections  1-3 
of  Chapter  2,  and  B.  A.  Krasovitskiy  wrote  Section  2 of  Chapter  1, 
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— temperatures  of  thawed  and  frozen  regions  respectively 
— initial  temperature  of  frozen  ground 
— radius  of  the  outer  wall  of  the  borehole 

— distance  from  the  thawing  front  to  the  axis  of  the  borehole 
— coefficients  of  thermal  .diffusivity  of  thawed  and  frozen 
ground  respectively 

— coefficients  of  thermal  conductivity  of  thawed  and  frozen 
ground  respectively 

— coefficients  of  heat  transfer  from  gas,  oil,  and  flushing 
fluid  to  the  wall  of  the  borehole 
— thickness  of  the  cement  collar  of  the  borehole 
— coefficient  of  thermal  conductivity  of  cement  stone 
— temperature  of  gas,  oil  in  well 
— latent  heat  of  fusion  of  water 
— moisture  content  of  the  ground  by  volume 
— time 

— time  scale 

— time  of  beginning  of  thawing 
— thickness  of  stratum  of  permafrost 
— specific  heat  of  gas,  oil 
— density  of  gas,  oil 
— central  speed  of  gas,  oil  in  well 
— well  output  and  flushing  fluid  velocity 
— pressure  at  entrance  to  stratum  of  frozen  ground 
— pressure  at  exit  from  stratum  of  frozen  ground 
— thickness  of  wall  of  pipe 

— coefficient  of  thermal  conductivity  of  wall  of  pipe 
— Joule-Thomson  coefficient 
— heat  equivalent  of  work 
— temperature  of  flushing  fluid  in  pipe 
— temperature  of  flushing  fluid  in  annular  space 
— coefficient  of  heat  transfer  between  flushing  fluid 
in  drill  pipe  and  in  annular  space 
— specific  heat  of  flushing  fluid 
— diameter  of  drill  pipe 

— temperature  of  flushing  fluid  at  entrance  to  drill  pipe 
— coefficient  of  compressibility 

— kinematic  viscosity  and  coefficient  of  thermal  conductivity 
of  flushing  fluid 


CHAPTER  1 

REVIEW  OF  COMPLETED  INVESTIGATIONS 


The  description  of  the  processes  which  take  place  during  the 
movement  of  gas  and  oil  in  pipes  is  based  on  the  use  of  three 
fundamental  laws:  the  laws  of  conservation  of  mass,  momentum,  and 

energy.  Analysis  of  these  processes  is  carried  out  using  the  methods 
of  continuum  mechanics. 

From  the  point  of  view  of  continuum  mechanics,  a liquid  (oil) 
and  a gas  are  defined  as  compressible  media  in  which  there  may  exist 
only  negative  normal  stresses  (pressure),  and  also  tangential 
stresses  due  to  viscous  friction.  In  this  case  all  three  conservation 
laws  may  be  given  in  the  form  of  several  mathematical  laws  which 
are  generally  valid  for  describing  the  motion  of  any  continuous 
medium. 

The  form  of  these  laws  is  simplified  considerably  for  one- 
dimensional motions  in  pipes  having  constant  cross-sections. 

Hydraulics  is  concerned  with  analyzing  such  motions.  Hydraulics  is 
one  of  the  applied  subfields  of  hydrodynamics;  the  methods  of 
theoretical  investigation  of  the  problems  of  hydraulics  proper  are 
continually  approaching  the  corresponding  methods  of  hydromechanics 
(methods  for  solving  boundary-layer  problems,  and  theories  of 
turbulence  and  gas  dynamics). 

The  fundamental  results  of  liquid  and  gas  pipe  hydraulics  may 


be  found  in  the  second  volume  of  the  survey,  "Mechanics  in  the 
USSR  over  50  years"  [1];  therefore,  in  what  follows  we  shall  only 
consider  those  investigations  which  are  directly  connected  with 
the  thermal  regime  of  wells. 

Section  1.  Investigations  of  Heat  Transfer  Processes  During  the 
Movement  of  Oil  and  Gas  in  Wells 

In  spite  of  the  obvious  generality  of  the  processes  which  take 
place  during  the  movement  of  a.  liquid  or  gas  in  pipelines  and  wells, 
the  study  of  the  thermal  regime  of  the  latter  was  only  begun  in 
1950;  moreover,  the  first  works  almost  totally  ignored  the  achievements 
of  one-dimensional  hydrodynamics  [2]  and  gas  dynamics  [3,  *0  connected 
with  the  calculation  of  flows  in  pipes.  This  is  explained  by  various 
fac  ors,  including  the  characteristics  of  the  operation  of  the  wells 
anu  the  variety  of  goals  which  are  set  during  the  investigation  of 
heat  exchange  in  wells  and  pipelines.  The  first  factor  will  be 
considered  below;  as  for  the  second  factor,  its  influence  is  clear, 
if  only  because  of  the  fact  that  the  first  investigations  of  the 
thermal  regime  of  wells  were  due  to  the  requirements  of  geothermal 
studies  and  geophysical  exploration  [5,  6]. 

In  a series  of  subsequent  works  [7,  8,  9],  questions  of  heat 
exchange  between  wells  and  rock  were  studied  in  connection  with  an 
investigation  of  the  effect  of  the  circulation  of  flushing  fluid 
during  boring  on  the  heat  field  of  the  ground  around  the  well. 

, Work  [7]  deserves  special  mention,  since  this  work  for  the 
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first  time  precisely  formulated  the  quasistationary  problem  of 
heat  e ange  between  a flow  of  liquid  in  a pipe  and  rock.  Sub- 
sequently this  approach  to  the  investigation  of  the  thermal  regime 
of  wells  was  used  rather  widely;  in  the  domestic  literature  its 
author  is  considered  to  be  I.  A.  Charnyy,  who  independently  carried 
out  an  analogous  investigation  [9]  several  years  after  work  [7]. 

The  basic  idea  of  works  [7,  9]  consisted  in  the  assumption 
regarding  the  significant  difference  in  the  time  scales  (if  one 
uses  the  contemporary  terminology  of  [10])  of  heat  exchange  in 
pipes  (forced  convection)  and  in  rock  (thermal  conductivity). 

Under  this  physically  obvious  assumption,  the  substationary  derivative 
in  the  energy  equation  may  be  neglected,  and  the  coefficient  of 
heat  transfer  from  the  liquid  to  the  rock  is  considered  to  depend 
weakly  on  time.  The  difference  between  works  [7]  and  [9]  consists 
only  in  the  choice  of  the  form  of  the  dependence  of  this  coefficient, 
on  time;  in  the  first  work  this  dependence  is  determined  from  an 
exact  solution  of  the  corresponding  thermal  conductivity  problem, 
and  in  the  second  work,  from  an  approximate  solution,  obtained  using 
the  method  of  quasistationary  states.  Subsequent  comparison  with 
the  exact  solution  showed  that  the  results  of  work  [9],  generally 
speaking,  are  suitable  for  sufficiently  large  values  of  time  [11]. 

The  ideas  of  work  [7]  are  found  also  in  the  investigation  of 
H.  Ramey  [12],  who  examined  the  thermal  regime  of  oil  and  gas  wells 
in  the  quasistationary  approximation.  Considering  the  gas  to  be 
ideal,  and  the  oil  to  be  an  incompressible  fluid,  he  obtained  from 
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the  equation  of  energy  conservation  an  ordinary  differential  equation 
of  the  first  order  for  determining  the  temperature.  The  dependence 
of  the  heat  transfer  coefficient  on  time  was  determined  from  an 
approximate  solution  of  an  axisymmetric  problem  for  a heat  conduction 
equation  suitable  for  large  values  of  time. 

As  was  indicated  by  V.  N.  Petrov  [131,  the  formulas  obtained 
by  H.  Ramey  give  sufficiently  accurate  results  when  calculating 
the  temperature  in  water  and  oil  wells  for  the  case  where  they 
operate  for  a sufficiently  long  period  of  time.  As  for  gas  wells, 
by  virtue  of  the  assumptions  used  (the  ideal  gas  assumption)  Ramey's 
formulas  are  valid  for  wells  that  are  not  deep  (on  the  order  of 
1000m)  and  for  pressure  gradients  that  are  not  large. 

An  attempt  to  take  into  account  the  real  properties  of  a gas 
was  undertaken  for  the  first  time  by  S.  A.  Bobrovskiy  and  V.  I. 
Chernikin  [14],  who  investigated  the  stationary  nonisothermal  movement 
of  gas  in  a well  taking  into  account  the  Joule-Thomson  effect  for 
the  case  of  a constant  mass  flow.  The  temperature  of  the  rock  was 
assumed  to  be  linearly  dependent  on  the  depth,  and  the  coefficient 
of  external  heat  exchange  and  the  oriflcing  coefficient  were  assumed 
to  be  constant.  In  addition,  it  was  assumed  that  the  pressure  was 
a linear  function  of  the  depth  of  the  well. 

Later  (J.  P.  Korotayev  analyzed  in  detail  the  stationary 
nonisothermal  flow  of  a real  gas  in  a well  for  the  case  of  a rock 
temperature  which  depended  linearly  on  depth.  The  results  of  his 
investigations  {153  amount  essentially  to  the  following. 
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For  the  equation  of  state  of  a gas  that  was  proposed  by  K.  V. 
Pokrovskiy,  integration  of  the  system  of  ordinary  differential 
equations  describing  the  motion  of  a real  gas  in  a well  may  be 
carried  out  using  the  method  of  successive  approximations,  where 
the  zeroeth  approximation  corresponds  to  the  case  of  an  ideal  gas. 

The  general  solution  is  constructed  in  the  form  of  a series  in 
powers  of  a small  parameter.  By  means  of  comparison  with  the  exact 
numerical  solution  it  is  shown,  that  for  the  majority  of  practically 
interesting  cases  the  first  approximation  is  sufficient.  However, 
even  in  this  case  the  solution  is  obtained  in  the  form  of  cumbersome 
quadratures,  which  compelled  the  author  to  look  for  the  possibility 
of  further  simplification.  Unfortunately,  the  method  selected  by 
him  is  not  always  convincing.  For  example,  it  is  pointed  out  that 
a taking  into  account  of  the  real  properties  of  the  gas  has  a greater 
effect  on  the  character  of  the  temperature  curve  than  on  the  pressure 
gradient  between  the  bottom  and  mouth  of  the  well.  However,  following 
this  the  temperature  is  averaged  over  the  depth  of  the  well  in  order 
to  simplify  the  method  of  calculation  of  the  pressure. 

It  is  necessary  to  note  that  the  results  obtained  by  U.  F. 
Korotayev  may  be  easily  extended  to  the  case  where  one  takes  into 
account  the  quasistationary  heat  exchange  between  the  gas  and  the 
rock. 

Further  investigations  of  the  thermal  regime  of  operating  wells 
are  connected  with  the  solution  of  nonstationary  problems  o f heat 
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exchange  between  a flow  of  liquid  or  gas  and  rock.  Froblems  of 
this  sort  belong  to  the  class  of  conjugate  problems  of  heat  exchange 
[10]. 

In  connection  with  the  study  of  the  thermal  regime  of  operating 
wells,  the  first  problem  of  this  sort  was  examined  in  1957  by  L.  Lesem 
and  coauthors  [16],  who,  for  the  case  of  an  ideal  gas  and  with  the 
aid  of  the  Laplace  transform,  solved  a system  of  linear  differential 
equations  in  partial  derivatives  that  consisted  of  the  energy  equation 
for  the  gas  in  the  well  and  the  heat  conduction  equation  in  rock. 

In  the  latter  equation,  heat  conduction  along  the  borehole  was 
neglected.  The  conditions  of  conjugation  between  the  heat  flows 
in  the  well  and  in  the  rock  corresponded  to  the  case  of  the  absence 
of  heat-transfer  resistance.  The  solution  was  obtained  in  the  form 
of  complicated  integrals  of  Eessel  functions,  the  graphs  of  which 
are  given  in  [16], 

The  results  of  L.  Lesem  et  al.  were  used  in  work  [17],  which 
was  devoted  to  a theoretical  investigation  of  the  thermal  regime 
of  water-injection  wells.  Comparison  of  the  results  of  the  calcula- 


tions with  the  data  of  the  conducted  experiments  proved  to  be 
satisfactory.  The  authors  of  work  [17]  note  the  difficulties  in 
calculating  the  quadratures  obtained  by  L.  Lesem  et  al.  in  [16], 
and  derive  simple  asymptotic  formulas  which  are  suitable  for  large 
values  of  time.  Later,  N.  A.  Avdonin  and  A.  A.  Buykis  obtained 
a formula  for  small  values  of  time  and  showed  that  it  corresponds 
to  the  solution  of  the  problem  of  the  injection  of  fluid  through 
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a gallery  [18].  Earlier  this  same  result  was  obtained  by  E.  E. 
Chekalyuk  using  a somewhat  different  approach  [19]. 

With  the  advent  of  high-speed  computers  there  arose  the 
possibility  of  a qualitatively  new  approach  to  the  investigation 
of  the  thermal  regime  of  wells.  This  possibility  was  successfully 
realized  by  V.  N.  Petrov  in  his  doctoral  dissertation,  completed 
in  1968.  Unfortunately,  the  fundamental  results  of  this  large 
and  interesting  work  were  published  in  a variety  of  low-circulation 
collections,  for  example  [20],  and  this  work  turned  out  to  be  less 
well-known  than  it  deserves  to  be. 

The  work  of  V.  N.  Petrov  is  primarily  concerned  with  an  analysis 
of  the  various  effects  caused  by  taking  into  account  the  real 
properties  of  a gas  and  by  the  effect  of  these  properties  on  the 
thermophysical  characteristics  of  the  gas.  The  calculations  carried 
out  by  him  showed  that  the  differences  in  the  mouth  temperature 
values  for  calculation  variants  distinguished  from  each  other  by 
the  equation  of  state  and  by  the  nature  of  the  dependence  of  the 
heat  capacity  on  pressure  and  temperature  attain  10-15°C  and  increase 
with  time.  All  the  calculations  were  made  for  a constant  mass  flow 
of  the  gas  and  under  several  physically  well-founded  assumptions. 

The  work  contains  a series  of  conclusions  that  are  important  from 
a practical  standpoint;  of  these,  we  note  the  results  relating  to 
the  integration  of  the  gradient-thermogram  of  gas  wells.  V.  N. 

Fetrov  showed  that  in  exploited  wells  the  gradient -thermograms  are 
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relatively  stable,  and  the  amplitudes  of  anomalies  attain  a maximum 
in  the  first  few  hours  after  the  well  is  started  up. 

In  recent  years  the  thermal  regime  of  gas  wells  has  attracted 
the  attention  of  investigators  in  connection  with  the  problem  of 
preventing  hydrate  formation.  This  problem  is  of  special  interest 
for  regions  where  long-terrfi  frost  occurs,  since  here  the  problem 
of  determining  the  thermal  regime  of  wells  becomes  especially 
difficult  due  to  the  presence  of  a moving  interface  between  thawed 
and  frozen  rock.  The  determination  of  the  equation  of  motion  of 
this  boundary  and  the  temperature  fields  in  the  thawed  and  frozen 
zones  in  itself  presents  considerable  mathematical  difficulties 
[21],  which  are  compounded  upon  the  simultaneous  examination  of  the 
heat  fields  in  the  rock  and  the  heat  flow  in  the  well.  It  is  not 
surprising  that  the  number  of  works  devoted  to  this  subject  is  very 
small . 

First,  let  us  mention  the  investigations  in  which  the  effect 
of  long-term  frost  is  taken  into  account  by  means  of  partitioning 
the  region  of  integration  into  two  zones  with  different  geothermal 
curves,  corresponding  to  ordinary  and  long-frozen  rock  [15,  22,  23]. 
It  is  obvious  that  such  an  approach  cannot  give  satisfactory  results, 
since  here  the  role  of  the  ice-to-water  phase  transition  is  not 
taken  into  account. 

B.  L.  Krivoshein  and  A.  A.  Koshelev  [24]  have  proposed  an 
algorithm  for  solving  the  problem  of  the  determination  of  the  tem- 
perature of  operating  wells  in  a region  of  long-term  frost  and  have 


announced  the  possibility  of  its  realization  on  a computer,  where 
it  is  implicitly  assumed  that  it  is  possible  to  separate  the  problem 
of  the  motion  of  ttie  gas  from  Stefan's  problem.  Works  [25,  26] 
formulate  and  solve  for  the  first  time  the  problem  of  the  heat 
exchange  between  an  operating  well  and  long-frozen  rock.  The  content 
of  these  works  will  be  presented  in  the  appropriate  section  of  this 
book . 

Section  2.  Methods  for  Solving  Problems  of  Heat  Exchange  Between 
Wells  and  Frozen  Rock 

The  intensity  of  the  heat  transfer  from  the  flow  in  a well 
into  the  ground  is  determined  by  the  nature  of  the  temperature  field 
of  the  ground  in  the  region  of  the  borehole.  During  heat  exchange 
between  the  flow  (which  has,  as  a rule,  a positive  temperature) 
and  the  surrounding  frozen  ground,  there  occur  in  the  latter  phase 
transitions  accompanied  by  the  absorption  of  a significant  amount 
of  heat.  This  leads  to  the  fact  that  the  rate  of  movement  of  the 
zero  isotherm  in  such  ground  sharply  decreases  in  comparison  with 
ground  in  which  there  are  no  phase  transitions,  and  the  intensity 
of  heat  transfer  for  the  corresponding  moments  in  time  increases. 
Thus,  in  order  to  calculate  the  temperature  of  a flow  moving  in  a 
well,  it  is  necessary  to  obtain  a solution  of  the  problem  of  the 
thawing  of  frozen  rock  in  the  space  around  the  well.  Very  many 
works  are  devoted  to  solving  this  problem.  In  this  paragraph  we 
shall  consider  only  the  most  significant  of  these  works. 

At  the  present  time,  the  majority  of  geocryologist s accept 
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the  following  definition  of  frozen  rock:  "Py  frozen  rock,  ground, 

and  soil  is  meant  rock,  ground,  and  soil  which  possess  a negative 
or  zero  temperature  and  in  which  at  least  part  of  the  water  has 
made  the  transition  into  the  crystalline  state  [27]." 

Upon  the  supply  of  a sufficient  amount  of  heat  to  the  frozen 
rock,  the  ice  cementing  the  separate  particles  of  rock  fully  or 
partially  turns  into  water.  This  sharply  worsens  the  cohesiveness 
of  the  rock,  which  then  serves  as  a reason  for  the  formation  of 
subterranean  cavities,  griffons,  gas  leakage  in  the  space  around 
the  well,  and  other  complications. 

The  water  in  the  frozen  rock  interacts  with  the  active  surfaces 
of  the  mineral  skeleton.  This  lowers  the  temperature  at  which  the 
water  freezes.  Simplifying  somewhat  the  rather  complex  physical 
picture  of  this  phenomenon,  one  may  say  that  the  thinner  the  capillary 
tube  in  which  the  water  is  located,  the  lower  its  freezing  tem- 
perature. Thus,  the  water  in  the  frozen  ground  freezes  in  some  | 

spectrum  of  temoeratures . Moreover,  it  may  be  confirmed  that  at 

i| 

any  negative  temperature  some  amount  of  the  water  in  the  frozen 
ground  will  be  in  the  unfrozen  state. 

On  this  basis,  A.  G.  Kolesnikov  [28]  proposed  the  following 
formulation  of  the  problem  of  the  freezing  of  the  ground.  The 
iciness  of  the  ground  i is  assumed  to  be  a known  function  of  the 
temperature  © , L.e.  i = i((pP).  The  thermophysical  properties  . | j 
of  the  ground  depend  on  the  iciness:  c = c(i);  ^ = ^(i).  A 

change  in  the  iciness  A leads  to  the  liberation  (absorption)  of 


the  latent  heat  of  fusion  ^aio  Au  , where  w is  the  moisture 
content . 

Thus,  the  heat  conduction  equation  for  freezing  moist  ground 
(for  the  one-dimensional  case)  assumes  the  form 
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(1.1) 


In  this  formulation,  the  heat  of  fusion  is  taken  into  account 
in  the  form  of  internal  sources  of  heat  that  are  distributed  over 
volume;  the  power  of  these  sources  depends  on  the  temperature. 

S.  H.  Cho  and  J.  E.  Sunderland  [29]  analogously  examine  the  problem 
of  freezing  or  thawing.  The  complexity  of  the  obtained  equations 
was  the  reason  that  the  solutions  of  the  problems  of  the  freezing 
or  thawing  of  the  ground  in  this  formulation  have  been  obtained 
only  for  the  simplest  cases  [30],  [31]-  Another  solution  of  the 
problem  of  freezing  or  thawing,  known  as  Stefan's  problem,  has  found 
much  wider  application.  " 

Study  of  the  iciness  curves  of  frozen  ground  shows  that  the 
overwhelming  part  of  the  water  in  the  ground  freezes  in  a spectrum 
of  negative  temperatures,  localized  around  0°C.  This  spectrum  is 
very  narrow  for  coarse  rock  and  is  somewhat  wider  for  fine  rock. 
Proceeding  on  the  basis  of  this  fact,  the  assumption  is  made  that 
all  the  phase  transitions  of  the  water  in  the  ground  proceed  at  0°C 
and  that  the  latent  heat  of  fusion  is  liberated  only  at  the  interface 
between  the  solid  and  liquid  phases.  In  this  formulation  the  problem 


of  the  thawing  of  frozen  ground  reduces  to  the  solution  of  the 
following  equations  (for  the  one-dimensional  case): 


The  problem  of  the  freezing  of  moist  ground  is  posed  analogously. 

The  first  solution  of  Stefan’s  problem  was  found  in  1831,  when 
G.  Lame  and  E.  P.  Clapeiron  [32]  found  a similitude  solution  for 
the  case  © „ = 0 . A similitude  solution  for  ® i-  0 was  given 

r & 

by  F.  Neumann  [33],  J.  Stefan  [3*0,  and  C.  Schwartz  [35].  N.  N. 

Verigin  [36]  found  a similitude  solution  for  the  axisymmetric 
Stefan's  problem  for  a linear  cold  source  (for  the  case  where  the 
ground  freezes).  Exact  solutions  of  separate  problems  were  obtained 
in  the  works  of  I.  G.  Fortnov  [37]  and  G.  A.  Tirskiy  [38].  These 
works  practically  exhaust  the  number  of  exact  solutions  of  problems 
of  the  Stefan  problem  type.  In  the  overwhelming  majority  of  cases 
(nonuniform  problems,  problems  with  cylindrical  symmetry,  problems 
with  various  boundary  and  initial  conditions,  etc.),  numerical, 
analog,  and  approximation  methods  are  used  to  solve  these  problems. 

Let  us  mention  the  most  important  of  these  methods. 

In  the  works  of  L.  I.  Rubinshteyn  [21,  39],  A.  Friedman  [40], 

J.  J.  Kolodner  [41],  and  other  authors,  there  was  developed  a method 
for  reducing  Stefan's  problem  tc  a system  of  integral  equations  of 
the  Volterra  type.  Solution  of  the  latter  is  carried  out  with  the 
aid  of  iterations  on  a computer. 

V.  G.  Melamed  [42]  developed  a method  for  reducing  Stefan's 
problem  to  a system  of  an  infinite  number  of  ordinary  differential 
equations.  Truncating  this  system  to  n equations  gives  an  approximate' 
solution,  which  as  n — * t»-0  approaches  the  exact  solution.  An 
advantage  of  the  method  is  the  possibility  of  using  standard  programs 
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for  integrating  the  obtained  system  on  a computer. 


A large  number  of  works  are  devoted  to  the  development  of 
finite-difference  methods  for  solving  Stefan's  problem.  As  an 
example,  we  may  mention  the  works  of  I.  Douglas  [43],  B.  M.  Budak 
[44],  and  G.  M.  Dussinberre  [45,  46]. 

Analog  methods  with  the  use  of  hydrointegrators  and  electronic 
differential  analyzers  have  found  wide  application  in  the  solution 
of  Stefan's  problem.  I.  M.  Kutasov  [47],  G.  I.  Man'kovskiy  [48], 

F.  Ya.  Novikov  [49],  and  others  investigated  the  temperature  fields 
around  freezing  and  heating  wells  using  hydrointegrators.  Hydro- 
integrators of  the  system  of  V.  S.  Luk'yanov  are  often  used. 

A study  of  the  formation  of  ice  with  the  aid  of  analog  computers 
was  made  by  R.  Howe  [50]  and  A.  L.  London  [51].  A problem  taking 
into  account  variable  thermophysical  properties  was  investigated  on 
an  electronic  differential  analyzer  by  D.  R.  Otis  [52], 

In  view  of  the  fact  that  an  exact  solution  of  Stefan's  problem 
is  connected  with  great  computational  difficulties,  various  approx- 
imation methods  have  found  widespread  application.  Let  us  examine 
the  most  important  of  these  methods. 

1.  The  method  of  the  successive  replacement  of  stationary 
states.  The  idea  of  the  method  belongs  to  L.  S‘.  Leybenzon,  who 
used  the  method  in  solving  a problem  concerning  the  hardening  of 
oil  in  a pipe  [53].  He  assumed  that  because  of  the  small  rate 
of  advancement  of  the  hardening  front,  the  temperature  distribution 
in  the  solid- phase  at  each  moment  in  time  differs  little  from  the 
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stationary  distribution  for  the  given  position  of  the  hardening 
front.  The  temperature  of  the  liquid  phase  was  assumed  to  be  equal 
to  the  pour  point. 


was 


Thus,  the  temperature  and  heat  flow  at  each  point  of  the 

solid  phase  were  determined  from  a solution  of  the  stationary 

problem  as  a function  of  the  coordinate  of  the  point  and  the 

coordinate  of  the  boundary.  The  obtained  expression  for  __ 

c * S 

substituted  into  Stefan's  condition,  which  led  to  an  ordinary 
differential  equation  for  S(t).  This  method,  because  of  its 
simplicity,  has  found  wide  application  in  engineering  calculations 
[54]. 

One  may  consider  the  work  of  F.  Kreith  and  F.  E.  Rorn.ie  [55] 
to  be  an  evolvement  of  this  method.  The  essence  of  this  work 
reduces  to  the  following.  Actually,  in  the  method  of  L.  S.  Leybenzon 
the  temperature  distribution  is  determined  from  equation  (1.2)  under 
the  assumption  that  the  inertia  term  is  ecluaI  to  zero. 

■ — ■■  (c ) 

In  [55]  this  solution  is  taken  as  the  zeroeth  approximation  (Q 
One  looks  for  a solution  to  the  problem  having  the  form 


- 22  e{'°. 


where 


is  determined  from  the  equation 
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The  problems  of  the  hardening  of  fluid  in  a pipe  and  in  a sphere 
for  a zero  initial  temperature  are  solved  in  [55]  using  this  method. 
The  method  is  analogous  to  the  one  which  was  used  by  M.  Ye.  Shvets 
[56]  and  Ye.  M.  Dobryshman  [57]  in  the  solution  of  the  equations 
of  a boundary  layer. 

2.  The  integral  method  in  its  contemporary  formulation  was 
developed  by  T.  Goodman  [58]  and  I.  Green  [59],  but  it  was  used 
earlier  by  I.  A.  Charnyy  [60]  in  solving  plane  and  axisymmetric 
Stefan  problems. 

The  essence  of  the  method  consists  in  replacing  the  heat  flew 
equation  by  the  heat  balance  relation.  The  latter  is  obtained 
by  equating  the  increment  in  the  heat  accumulated  by  the  body  per 
unit  time  to  the  heat  flow  through  the  moving  boundary.  T.  Goodman 
obtains  the  relation  formally,  by  Integrating  the  heat  flow  equation 
over  the  spatial  coordinate  within  the  limits  of  a single  phase, 
and  I.  A.  Charnyy  obtains  it  from  physical  considerations. 

Next,  the  temperature  profiles  are  presented  in  the  form  of 
polynomials  with  coefficients  depending  on  the  coordinate  of  the 
boundary  (for  the  linear  problem),  or  in  the  form  of  a logarithmic 
profile  (for  the  axisymmetric  problem).  In  order  to  obtain  such 
a profile  in  a semi-infinite  region,  one  introduces  the  concept  of 
the  zone  of  thermal  influence,  at  the  boundary  of  which  the  heat 
flow  is  assumed  to  be  equal  to  zero,  and  the  temperature,  equal  to 
the  unperturbed  temperature  of  the  body.  The  obtained  temperature 
profiles  are  substituted  into  the  heat  balance  relation  and  into 


Stefan's  condition,  which  leads  to  a system  of  ordinary  differential 
equations  with  respect  to  the  unknown  functions. 

Kh.  R.  Khakimov  [6l],  using  this  method  to  solve  the  problem 
of  the  freezing  of  the  ground  around  a freezing  borehole,  assumed, 
in  order  to  simplify  the  problem,  that  the  radius  of  the  front  of 
thermal  influence  was  proportional  to  the  radius  of  the  thawing 
front.  The  coefficient  of  proportionality  was  determined  by  him 
experimentally.  V.  A.  Maksimov  et  al.  [62]  used  the  integral  method 
in  solving  the  problem  of  the  melting  of  a semi-infinite  red  in  a 
medium  having  a constant  or  varying  temperature,  with  ablation  or 
without  it.  In  the  aforementioned  work  [29]  the  integral  method 
was  used  in  solving  a problem  with  an  extended  phase  transition 
front.  T.  Goodman  [63]  used  this  method  in  solving  Stefan's  problem 
for  a body  with  variable  thermophysical  properties. 

3.  The  method  of  continuations.  A significant  feature  of  this 
method  is  the  fact  that  a fictitious  body  is  introduced,  the  shape 
of  which  is  invariant  with  time  and  coincides  with  the  shape  of 
the  body  in  question  up  to  the  beginning  of  the  phase  transformations 
Thus,  in  connection  with  problem  (1.2)-(1.7),  the  region  in  question 
is  0 x d 00  , in  which  equation  (1.3)  is  satisfied. 

In  order  that  the  conditions  be  satisfied  at  the  moving  boundary, 
a fictitious  heat  flow  is  introduced  at  the  fixed  boundary.  After 
satisfying  the  boundary  conditions,  one  arrives  at  an  integrodif f eren 
tial  equation  for  determining  the  heating  capacity  of  the  fictitious 
heat  flow.  This  equation  is  solved  numerically  or  by  expanding  the 
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desired  solution  in  a series.  Such  a method  was  developed  by  E.  A. 
Boley,  who  used  this  method  in  solving  one-dimensional  [64]  and 
axisymmetric  [65]  Stefan  problems  with  ablation  and  without  it. 

V.  I.  Antipov  and  V.  V.  Lebedev  [66],  having  solved  problem 
(1.2)— (1.7)  with  variable  initial  and  boundary  conditions,  introduced 
two  fictitious  temperatures:  one,  depending  on  time,  into  condition 

(1.4);  the  other,  depending  on  the  coordinate,  into  condition  (1.5). 
After  satisfying  the  conditions  at  the  moving  boundary,  they  arrived 
at  a system  of  two  integrodif f erential  equations,  which  were  solved 
by  expanding  the  unknown  functions  in  a series. 

G.  A.  Martynov  [67]  used  this  method  in  solving  the  inverse 
Stefan  problem. 

4.  The  method  of  contour  integrals.  The  method  was  proposed 
by  G.  A.  Grinberg  [68]  for  solving  problem  (1.2)  — (1.7)  for  (S)p  = 0. 
The  new  variable  z = x - S(F)  is  introduced.  Then  equation  (1.2) 
acquires  the  form 

4.  . A = (i.8) 

< II  0?  ' 


The  function 


exp[W-.|I(z+S(£-))]i 


where  'X  is  any  complex  constant,  is  a particular  solution  of 
equation  (1.8)..  The  general  solution  may  be  given  in  the 


form 


wi  = -TT  i c (>-)  exp  f >J—  I >.  i r + 5 (0 ']  i//., 
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where  c ( ^.  ) is  an  arbitrary  function.  The  contour  of  integration 
P is  chosen  in  the  form  of  a line  parallel  to  the  imaginary  axis 
in  such  a way  that  all  the  singular  points  of  the  integrand  remain 
to  the  left  of  the  line.  From  condition  (1.6)  we  obtain 

O --  j 

i c<?-)  e*P  [>./'—  I a,  —3  (?)]  d\  = 0.  (1.10) 
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From  this  relation  the  form  of  the  function  c ( \ ) is  determined. 

The  method  is  used  advantageously  when  solving  the  inverse 
Stefan  problem,  i.e.  when  determining  the  temperature  fields  from 
a given  equation  describing  the  movements  of  the  boundary  [67,  69], 
but  this  method  may  also  be  used  for  solving  the  direct  problem 
under  the  conditions  of  a constant  initial  temperature  and  a constant 
temperature  at  the  boundary. 

5.  The  variational  method  (Biot's  method).  M.  A.  Biot  is 
justly  considered  to  be  the  author  of  this  method,  in  spite  of  the 
fact  that  the  variational  formulation  of  the  theory  of  heat  conduction 
was  first  proposed  by  Rosen  [70].  The  point  is  that  M.  A.  Biot  not 
only  developed  the  technique  for  applying  the  variational  principle 
to  dissipative  processes,  but  also  derived  general  differential 
equations  of  the  Lagrange  type  for  a thermodynamic  system  subjected 
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to  irreversible  changes  [71].  Later,  L.  G.  Chambers  used  this 
technique  for  solving  heat  conduction  problems  [72]. 

The  essence  of  the  method  of  M.  A.  Biot  consists  in  the 
following  [73]-  One  introduces  the  heat  flow  vector  field 

Q = Q(v„  x,). 

where  Q is  the  amount  of  heat  per  unit  area  that  has  passed  through 
the  given  cross  section  in  a given  direction  since  the  moment  the 
process  began;  s^  are  the  generalised  coordinates;  and  x.,  are  the 
spatial  coordinates.  The  following  invariants  are  introduced: 
the  thermal  potential 

. u = -4-  f cii&'dv 


and  the  dissipative  function 


Here 
Q = 


Then 


v is 

?a 

diz 

the 


the  volume 
; and  \ 
variat ional 


of  the  body  in  question;  Q is  the  temperature 
is  the  coefficient  of  thermal  conductivity, 
principle  is  formulated  in  the  following  manner 


bu  + 60  = j QbQiir, 


(1.11) 


where  F is  the  surface  of  the  body.  This  relation  is  equivalent  to 
the  heat  conduction  equation. 

If  one  introduces  the  concept  of  the  thermal  force 

vuu5 s,  - \MQJF, 
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then  the  variational  principle  may  be  written  in  the  form  of 
Lagrange's  equation: 


In  order  to  apply  this  principle  to  Stefan's  problem,  the 
latent  heat  of  fusion  is  introduced  in  place  of  Q,  and  the  coordinate 
of  the  moving  boundary  is  used  as  the  generalised  coordinate. 

Given  a temperature  profile  in  the  form  of  a polynomial,  one  arrives 
at  an  ordinary  differential  equation  for  S(F) . 

As  is  apparent  from  what  has  been  presented,  there  exist  a 
vast  number  of  methods  for  solving  Stefan's  problem;  these  methods 
have  been  developed  mainly  for  linear  and  axisymmetric  cases.  For 
the  solution  of  engineering  problems  concerning  the  thawing  of 
frozen  ground  around  a well  in  a plane  fornulation,  the  method  of 
successive  approximations  is  most  acceptable.  This  method  combines, 
as  will  be  shewn  below,  simplicity  and  a high  accuracy  of  the  obtained, 
solutions.  At  the  same  time,  practically  no  approximation  methods 
exist  for  solving  two-  and  three-dimensional  Stefan  problems,  which 


arise,  in  particular,  when  calculating  the  thermal  regimes  of  wells 
bored  and  exploited  in  a region  of  long-term  frost. 


CHAPTER  2 

HEAT  EXCHANGE  BETWEEN  FLOWS  IN  A WELL  AND  ROCK 


Section  1.  Derivation  of  a System  of  Equations  for  the  Nonisothermal 
Motion  of  Liquid  and  Gas  in  Pipes 

In  accordance  with  the  fundamental  assumptions  of  continuum 
mechanics,  for  the  description  of  the  nonisothermal  flow  of  a 
compressible  fluid  in  a well  we  use  three  conservation  laws: 
conservation  of  mass,  momentum,  and  energy. 

1.  The  law  of  conservation  of  mass  or  the  equation  of  continuity 
in  the  case  of  the  presence  of  internal  sources  of  mass  has  the 
form  [7^] 


+ div  (pi')  = m, 

0 1 


(2.1) 


where  is  the  density;  V is  the  velocity  vector;  m is  the  rate 

of  increase  in  the  mass  per  unit  volume  due  to  internal  sources. 

In  what  follows  we  shall  everywhere,  with  the  exception  of  special 
cases,  consider  that  m = 0. 

In  the  case  of  one-dimensional  motion  in  pipes  having  constant 
cross  sections,  from  (2.1)  we  obtain 
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In  practice  wells  most  often  operate  with  a constant  output. 


which  corresponds  to 


fvrf  - M 


c 


and 


(2.3) 


(2.4) 


Here  f is  the  cross  section  of  the  pipe,  and  M is  the  mass  flow. 

For  an  incompressible  fluid  (^  = const),  from  (2.2)  it  follows  that 

* 0?  (2.5) 


2jt 

9* 


or 


V * v i-t)  . (2,6) 

2.  The  dynamic  equation.  In  order  to  derive  the  equation  the 
following  theorem  is  used:  the  change  in  the  main  vector  of  the 

momenta  of  a system  of  material  points  is  equal  to  the  impulse  of 
all  the  external  mass  and  surface  forces  [7^1- 

The  application  of  this  theorem  in  pipe  hydraulics  gives  [75] 

— + oT  ^ f °')  = - Pg  sin  a - — p |?| v,  (2.7) 
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where  p is  the  pressure;  g is  the  acceleration  due  to  gravity;  &C  is 
the  angle  of  inclination  of  the  pipe,  measured  from  the  horizontal; 

\ is  the  coefficient  of  hydraulic  resistance;  and  D is  the  diameter 
of  the  pipe. 

For  wells  one  may  take  sir.  oc  = 1.  If  one  uses  (2.3),  then 
from  (2.7)  we  obtain 


Here  z(x)  is  the  ordinate  of  the  axis  of  the  pipe,  measured  from 
the  horizontal  plane;  I is  the  mechanical  equivalent  of  heat;  u is 
the  unit  internal  energy;  q is  the  heat  flow  through  the  wall  of 
the  pipe;  and  k is  the  coefficient  of  thermal  conductivity  of  the 
gas . 

Let  us  transform  the  left  hand  side  of  equation  (2.10)  in 
the  following  manner: 


P (gz  + E\ ] + ^[py [gz + + £jj  = p 4r(Sz  + £J  + 

+ (s* + E)  f +-P'J  3T  (s*  +E  + f)  + (s'z  + E + fy-^r 

-p£(*;  + £)  + P*£(f)  + f 


(r. 


-/«  + 


Now  equation  (2.10)  may  be  put  into  canonical  form  [74]: 


- Pgz'v  - 


0 fpv) 
ox 


+pQ . 


(2.11) 


where 


In  the  majority  of  practically  interesting  cases  the  heat 
conduction  may  be  neglected,  and  then  the  equation  of  energy  during 
motion  in  the  well  acquires  the  form  [next  page]: 


P jp Pgv- 


(2.12) 


(p~) 

Ox 


In  the  general  case  the  system  of  equations  (2.2),  (2.8),  and 
(2.12),  augmented  by  the  equation  of  state 

p * ^ (2-13) 

is  a system  of  equations  of  the  hyperbolic  type,  possessing  three 
real  characteristics  [75]: 


tlx 

at 


v. 


dx' 

(ft 


V ± c, 


(2.1*4) 


where 


is  the  speed  of  sound  in  the  gas. 

If  the  external  heat  flow  is  a given  function  of  x and  t,  then, 
as  is  shown  in  [75],  system  (2.2),  (2.8),  (2.12),  and  (2.13)  may  be 
reduced  either  to  characteristic  form  and  integrated  by  the  method 
of  characteristics,  or  to  a system  of  finite-difference  equations 
based  on  an  implicit  absolutely  stable  difference  scheme  [77]. 

The  matter  is  more  complicated  in  the  case  where  the  external 
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heat  flow  depends  on  the  nature  of  the  heat  exchange  between  the 
flow  in  the  well  and  the  rock.  Then  to  the  given  system  of  equations 
it  is  necessary  to  add  an  equation  describing  the  propagation  of 
heat  in  the  rock,  and  conditions  of  conjugation  of  the  heat  flows. 

For  rock  having  constant  thermophysical  properties  this  equation 
will  be  Fourier's  equation 


It  is  natural  to  assume  that  the  heat  from  the  well  propagates 
mainly  in  the  radial  direction,  which  allows  one  to  write  the  last 
equation  in  the  form 

_je  /()-h  _i_  de\ 

c n 1 ( d'r-  + r • $ )'  (2.15) 

where  Q)  is  the  temperature  of  the  rock,  and  a^  is  the  coefficient 
of  thermal  diffusivity.  At  the  contact  between  the  rock  and  the  well 
one  must  satisfy  the  boundary  condition 

for  r = , (2.16) 

where  \ ( is  the  thermal  conductivity  of  the  rock;  is  the 

total  heat  transfer  coefficient;  and  a is  the  radius  of  the  borehole. 
In  what  follows  it  is  convenient  to  operate  with  the  temperature  of 
the  gas;  therefore,  let  us  transform  the  energy  equation  (2.12). 
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Reading  off.  from  (2.12)  the  equation 


-pgu 


0 (PV) 
Ox 


+ N in, 


(2.17) 


expressing  the  theorem  regarding  the  change  in  the  kinetic  energy 
in  differential  form  [7^],  we  obtain 


o 


(<1— Min, 


(2.18) 


where  is  the  power  of  the  internal  forces  per  unit  volume. 

For  the  case  of  one-dimensional  motion  of  the  gas. 


Mi, 


Ap * 

r Ox  * 


which,  taking  into  account  the  equation  of  continuity,  may  be  put 
in  the  form 


Mi„  = - 


Ap  dp 
P dl  ' 


(2.19) 


If  we  now  pass  from  the  internal  energy  to  the  enthalpy 


(2.20) 


and  express  q In  terms  of  the  heat  flow  at  the  wall  of  the  borehole 


(2.21) 
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then  equation  (2.l8)  acquires  the  form 


dt  it  dr  La 


(2.22) 


Using  the  expression  for  the  enthalpy  of  a gas 


(ii=cf(dr—n.dp), 


(2.2?) 


where 


A 


is  the  Joule-Thomson  coefficient,  we  obtain 


<if  , , <.'/>  2 ae 

PC„-T  - (PCPH,  + A)-jr  = — A,  — r 


<u 


or,  passing  to  partial  derivatives  and  assuming  that  the  velocity 
and  pressure  of  the  gas  do  net  depend  on  time, 


ar  , ar  , , ap  2 ae 

pcP + pcpr  - «/ (pCpfi,-  + A)  ^ - — X,  -rr- 


(2.24) 


The  latter  assumption  is  based  on  the  fact  that  the  redistribution 
of  the  pressure  proceeds  much  more  rapidly  than  the  redistribution 
of  the  temperature,  and  consequently,  if  the  well  operates  with  a 
constant  output,  then  a very  short  time  after  the  well  is  started 
up  it  enters  into  the  quasistationary  regime,  in  which  the  slow 
changes  in  the  pressure  and  velocity  over  time  will  be  due  to  heat 
exchange  with  the  rock  and  to  the  effect  of  orificing  [20]. 

If  now  we  introduce  the  mass  flow  of  the  gas  in  accordance  with 
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dc 


24) 

we 

obtain 

( Ll : + 

A 

■')  cpM  = 2rtai1  ~ 

, 

V 
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j P Ox  1 Or 

\r  =a 

(2.25) 


For  the  case  of  a constant  mass  flow,  the  system  of  equations 
(2.9),  (2.13),  and  (2.25)  must  be  augmented  by  two  boundary  conditions 
and  one  initial  condition.  Usually,  the  pressure  and  temperature 
at  the  bottom  are  given  as  boundary  conditions.  At  the  same  time, 
it  is  clear  that  these  values  must  be  determined  from  the  solution 
of  the  corresponding  problem  of  the  nonisothermal  filtration  of  the 
gas  or  using  the  well-known  relations  of  subterranean  hydrodynamics, 
if  one  is  talking  about  a liquid  (oil). 


Section  2.  Determination  of  the  Temperature  of  Gas  at  the  Bottom 
of  Wells 

In  order  to  obtain  a complete  description  of  the  processes 
which  occur  during  the  filtration  of  a liquid  or  a gas  in  a porous 
medium,  it  is  necessary,  analogous  to  what  was  done  above,  to  derive 
a system  of  equations  of  continuity,  motion,  and  energy.  Such  a 
derivation  is  carried  out  in  a series  of  works  [19,  78-80].  However, 
the  obtained  system  of  equations  is  so  complex  that  the  solution  of 
the  corresponding  boundary  problems  turns  out  to  be  possible  only 
after  a series  of  simplifications  or  with  the  use  of  numerical  methods 
This  is  not  surprising,  if  one  recalls  that  the  solution  of  the 
problems  of  isothermal  filtration  in  itself  presents  considerable 
difficulties . 


**  ^ ^ 


For  these  reasons, 


analysis  of  the  temperature  fields  during 


r 


the  filtration  of  a liquid  or  gas  is  made  mainly  under  the  assumption 
that  the  effect  of  the  temperature  on  the  nature  of  the  redistribution 
of  the  pressure  is  negligibly  small  [19,  78].  With  such  an  approach 
it  is  considered  that  the  pressure  field  is  known  from  the  solution 
of  the  corresponding  problem  of  subterranean  hydrodynamics,  and 
once  it  is  so,  substitution  of  the  necessary  relations  into  the 
energy  equation  allows  one  to  find  in  principle  the  temperature  of 
the  gas  or  liquid. 

The  assumption  mentioned  above  is  based  on  the  fact  that  the 
rate  of  change  of  the  pressure  is  much  greater  than  the  rate  of 
change  of  the  temperature.  This  actually  takes  place  during  the 
filtration  of  an  elastic  liquid  in  strata  with  good  collector 
properties,  but  is  not  necessarily  true  during  the  filtration  of 
gas  in  deeply  located  strata  with  low  permeabilities.  For  example, 
for  stationary  filtration  of  an  ideal  gas,  from  Darcy's  law  and 
from  the  system  of  the  equations  of  continuity 

div(fa?)  =0  (2.26) 

and  energy 

_ ..  (227) 

div grad  T—cpyj  grad  F=0  ' ' 

it  follows  that  this  process  is  described  by  Laplace's  equation  for 


the  function 


u=pt+2bT, 


(2.28) 


where  p,  T are  the  dimensionless  pressure  and  temperature  respectively; 


C f 


r:  1 . 

* 1 l c j 

fi  T 


(2.29) 


'J'  is  the  bulk  density  of  the  gas  at  pressure  pc  and  temperature 
T0;  m is  the  porosity;  X,  is  the  coefficient  of  piezoconductivity 
for  an  initial  stratum  pressure  of  p^;  and  is  the  initial  stratum 
temperature . 

By  analogy  with  the  problems  of  thermoelasticity,  let  us  call 

b the  coupling  parameter.  In  the  given  case  it  characterizes  the 

Increase  in  pressure  due  to  coupling  of  the  velocity  and  temperature 

fields.  Estimates  show  that  under  the  usual  conditions  this  parameter 

-2  _■? 

is  on  the  order  of  10  to  10  . However,  for  strata  with  good 

thermophysical  and  poor  collector  properties,  it  may  have  a considerably 
higher  value. 

Thus,  for  the  case  of  a small  coupling  parameter,  for  determining 
the  temperature  at  the  bottom  of  boreholes  one  may  use  the  approximate 
formula  of  E.  B.  Chekalyuk  [19]: 


G * rb 


AT  = Hi- 


-\p  , , 

IT ln  V 


(2.30) 


where 


% ■ \ w 


is  the  radius  of  the  zone  of  influence  of  the  borehole;  r is  the 
radius  of  the  borehole;  h is  the  power  of  the  stratum;  c„  is  the 
volumetric  heat  capacity  of  the  rock;  G is  the  mass  yield  of  the 
well;  and  t is  the  time. 

If  the  dependence 


is  determined  using  the  well-known  formula  of  E.  E.  Chekalyuk  [19], 
then  after  straightforward  algebra  expression  (2.30)  acquires  the 
form 

^ ~ rb  A r = 

- cr 

Calculations  using  formulas  ( 2 . 30 )- ( 2 . 31 ) m ay  be  made  for  the 
case  of  small  depressions.  If  the  depressions  are  great,  then  this 
invalidates  the  assumptions  regarding  the  constancy  of  the  thermo- 
dynamic characteristics  of  the  system  that  were  made  in  deriving 
these  formulas. 

In  this  case,  let  the  coupling  parameter  again  be  small.  T 
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• the  pressure  may  be  determined  using  the  well-known  formulas  of 
subterranean  hydraulics  obtained  for  isothermal  filtration,  for 
example  using  the  formula  [79] 


(2.32) 


where  ^ is  the  viscosity  of  the  gas,  and  k is  the  permeability, 
and  then,  using  the  energy  equation  for  the  case  of  negligibly 
small  heat  conducting  flows 


•ir  , dp 

T + Hi  = 0, 
dr  dr 


(2.33) 


one  may  determine  the  temperature  from  the  solution  of  the  ordinary 
differential  equation  of  the  first  order 


7^  + [r-  p[r)]f  [r)  = 0, 


(2.3*0 


obtained  by  substituting  relation  (2.32)  into  equation  (2.33). 

If  the  coupling  parameter  cannot  be  considered  to  be  negligibly 
small,  then  it  is  necessary  to  determine  the  temperature  from  the 
solution  of  a system  of  ordinary  nonlinear  differential  equations, 
containing  besides  (2.33)  the  equation  of  continuity 


div(fti')  —0, 


(2.35) 
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the  equation  of  motion  (Darcy's  law) 


-fr  jp 


and  the  equation  of  state  of  a real  gas 

~ - z (, P , T)  RT. 


(2.36) 


(2.37) 


Combining  these  three  equations,  we  obtain 

^£-Z(p.T)T.  (2.38) 

The  system  (2.33),  (2.38)  is  solved  on  a computer  using  standard 
procedures.  The  fundamental  difficulty  connected  with  this  is  the 
calculation  of  the  functions  M?-  r)  and  l(  p>r)  . 

Attempts  to  select  two-dimensional  polynomials  using  tabulated  data 

have  so  far  failed  to  give  satisfactory  results  [81].  There  are 

also  problems  in  introducing  these  tables  into  the  memory  of  a computer 

[13]. 

It  is  probably  most  sensible  from  a practical  point  of  view  to 
use  various  semi-empirical  equations  of  state.  For  example,  in 
work  [75]  in  analyzing  the  nonisothermal  motion  of  a gas  in  a pipe, 
the  Berthelot  equation  was  used.  In  work  [15],  where  the  nonisothermal 
stationary  flow  of  gas  in  a well  was  studied,  an  equation  with 
three  virial  coefficients  was  used.  In  both  cases  the  calculations 
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gave  satisfactory  results. 

It  is  necessary,  however,  to  remember  that  the  use  of  any  of 

the  equations  of  state  is  valid  only  in  some  range  of  pressures 

and  temperatures;  moreover,  the  best  method  of  verification  is 

the  construction,  using  the  chosen  equation,  of  the  inversion  curve 

corresponding  to  the  case  where  the  Joule-Thcmson  coefficient  it 

equal  to  zero  [31].  In  particular,  such  a verification  has  shov/n 

that  the  applicability  of  the  Eerthelot  equation  is  restricted  to 

the  value  of  the  derived  temoerature  T,  =1.3,  and  an  equation 

der 

with  virial  coefficients  in  general  poorly • describes  the  inversion 
curve  [ 8 1 1 . 

Section  3.  The  Effect  of  the  Heat  Field  of  Rock  on  the  Thermal 
Regime  of  Wells 

Earlier  (Chapter  1 and  Section  1 of  Chapter  2)  it  was  mentioned 
that  the  effect  of  nonstationary  heat  exchange  between  a flow  in  a 
well  and  rock  may  be  taken  into  account  using  the  method  of  successive 
replacement  of  stationary  states  [9,  191.  A different  idea  was 
enunciated  by  E.  B.  Chekalyuk  [19].  He  assumed  in  deriving  an 
equation  describing  the  nonstationary  distribution  of  the  temperature 
in  a well  that  the  flow  of  heat  through  the  wall  of  the  pipe  may 
be  written  in  the  form 


l 
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where  c*.  (_tO  is  the  time-dependent  dimensionless  heat  exchange 
coefficient,  and 

AT  « T - 

where  the  rock  temperature  (5y  no  longer  depends  on  time. 

A large  number  of  works  are  known  where  various  special  cases 
of  formula  (2.39)  are  used.  In  the  general  case,  in  order  to  apply 
(2.39)  it  is  necessary  to  know  the  form  of  the  function  x fr)  ■ 

In  practice  the  formulas  of  I.  A.  Charnyy  [9,  79]  and  E.  B.  Chekalyuk 
[19]  are  most  often  used;  these  formulas  are  obtained  from  an 
approximate  solution  of  the  corresponding  heat  conduction  problems. 

Below  it  will  be  shewn  that  the  intuitively  introduced  coefficient 
<*(*)  has  a simple  physical  meaning.  It  is  simple'st  to  do  this  using 
the  example  of  the  operation  of  an  oil  well.  Let  the  oil  well 
operate  with  a constant  yield  Q.  The  temperature  of  the  rock  linearly 
increases  with  depth.  For  simplification  of  the  calculations  let 
us  assume  that  the  heat -transfer  resistance  of  the  walls  of  the 
borehole  is  negligibly  small.  Then  the  system  of  equations  describing 
the  given  process,  after  a series  of  algebraic  manipulations,  has 
the  form: 


■nr  , dsr  _ ^ or_.  (2.^0) 
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are  the  heat  capacity  and  density  of  the  rock; 
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(2.42) 


(2.43) 


(2.44) 


(2.45) 


) 

and  P 


the  geothermal  gradient. 


The  given  boundary  and  initial  conditions 


correspond  to  the  equality  of  the  temperatures  in  the  well  and  in 
the  rock  before  the  well  begins  operating,  and  also  to  the  equality 
of  the  stratum  and  bottom  temperatures  during  the  operation  of  the 
well . 

From  the  solution  of  equation  (2.41)  with  boundary  conditions 
(2.42)  and  (2.45)  and  initial  condition  (2.43)  and  using  the 
convolution  theorem  [82],  we  obtain 

r 

(tf)  , = -TT.f  T(r-i)Af(i)</e,  (2.46) 

\ //•=!  Q 


where 
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(2.47) 


JQ(u)  and  Yq(u)  are  Bessel  functions  of  the  first  and  second  kind, 
of  the  zeroeth  order.  Substituting  (2.47)  into  (2.40),  we  obtain 


o\t  ii_vr 
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or 


= 0. 


(2. 48) 


The  function  ' if)  (_7r)  has  a simple  physical  meaning.  It  is 
equal  to  the  dimensionless  heat  flow  at  the  wall  of  the  borehole 
[83].  Using  the  Laplace  transform  one  may  show  that  fcr  small  values 
of  time 

fp(-r)  « l/y-TT,  (2.49) 
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and  then  equation  (2.48)  coincides  with  the  equation  obtained 
by  E.  B.  Chekalyuk  for  the  flow  of  liquid  in  a plane  channel  [19]. 
This  confirms  the  correctness  of  the  result  obtained  here,  since 
it  is  known  [83]  that  for  small  values  of  time  the  solutions  of 
axisymmetric  and  plane-parallel  problems  coincide. 

Comparing  the  third  term  of  equation  (2.48)  and  formula  (2.3S), 
we  convince  ourselves  that  the  function  y'Ct)  corresponds  to  the 
function  (t)  , i.e.  this  is  also  the  dimensionless  heat  flow 
through  a cylindrical  surface  over  which  is  maintained  a constant 
temperature . 

In  precisely  the  same  way  one  may  shew  that  in  the  presence 
at  the  wall  of  the  borehole  of  a heat-transfer  resistance,  i.e. 
for  the  replacement  of  boundary  condition  (2.42)  by  a boundary 
condition  of  the  third  kind,  equation  (2.48)  retains  its  form, 
with  the  exception  of  the  fact  that  for  the  function  (f  £t)  it  is 
necessary  to  read  in  place  of  (2.47) 


1 r n{  ' u ! « ;|"^i  Vo  i- Nu  r l"' t (n)  N'u  yt  (i/)|sj  * 


(2.50) 


where  (u)  and  Y^(u)  are  Bessel  functions  of  the  first  and  second 
kind,  of  the  first  order,  and 


is  the  Nusselt  number. 
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It  is  me st  convenient  to  find  the  solution  of  the  various 
boundary-value  problems  for  equation  (2.48)  using  the  Laplace  transform 
As  an  example,  let  us  examine  the  solution  of  this  equation  tor 
boundary  condition  (2.44)  and  initial  condition  (2.43).  Using 
the  Laplace  transform 

cv 

«•  - | r)dT 

6 

in  place  of  (2.43),  (2.44),  and  (2.48),  we  obtain 
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(2.51) 
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(2.52) 


where  p is  the  transform  parameter; 


T (p)  = j c-^T  (T)dT. 


(2.53) 


The  solution  of  equation  (2.51)  for  boundary  condition  (2.52) 
has  the  form 


qi- 


|p  4-  P(MP)!  P 
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In  the  general  case,  in  order  to  pass  from  the  representation 
of  w to  the  original  quantity  T it  is  necessary  to  use  one  of 

the  numerical  methods  of  [82],  as  is  done  in  work  [11]  in  calculating 
the  thermal  regime  of  bored  wells.  The  solutions  for  small  may 

be  found  by  using  the  asymptotic  function  , and  for  large 

values  of  time  the  solutions  may  be  found  using  the  method  presented 
in  [833- 

For  many  practically  interesting  cases  the  method  of  quasi- 
stationary states  assures  sufficient  accuracy.  In  accordance  with 
this  method,  the  solution  of  the  just-examined  problem  will  have 
the  form 


ST  = [ 1 — g-CKTUl 


(2.55) 


where  according  to  [19] 


m (|  -j- 


(2.56) 


If  the  temperature  at  the  bottom  of  the  borehole  is  less  than 
the  stratum  temperature  due  to  orificing,  then  the  corresponding 
solution  has  the  form 


AT  = A7\c~a(t>*  + 


nu.a  (x) 


(2.57) 


From  (2.57)  one  may  find  the  depth  at  which  the  gas  that  is 
cooled  as  a result  of  orificing  is  heated  up  to  the  temperature 
of  the  rock  at  the  given  point.  Setting  the  left  hand  side  equal 
to  zero,  we  obtain  after  straightforward  algebra 


5 » (2.58) 

For  the  majority  of  practically  interesting  cases  the  second  term 
under  the  logarithm  sign  is  on  the  order  of  10-^.  Then  from  (2.58) 
it  follows  that 

= — I 'Q  ■ 

Taking  into  account  the  expression  for  , we  finally  obtain 

(2.59) 

where  x is  the  distance  from  the  bottom.  A taking  into  account  of 
the  real  properties  of  the  gas  during  its  motion  in  the  well  leads 
to  an  increase  in  this  distance. 

Section  Formulation  of  the  Problem  of  the  Heat  Exchange  Eetween 
Frozen  Rock  and  a Well 

Let  us  examine  the  thermodynamical  regime  of  exploitation  of 
a well  in  frozen  ground  (Fig.  1).  It  is  necessary  to  note  that 
strata  of  frozen  ground  practically  never  attain  the  designation  of 
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producing  beds,  having,  as  a rule,  a very  high  positive  temperature. 
Consequently,  the  cross  section  of  a well  in  a region  of  permafrost 
may  be  presented  in  the  following  manner:  the  producing  ted,  a strat 

of  warm  (unfrozen)  ground,  and  a stratum  of  frozen  ground.  During 
the  motion  of  gas  (oil)  in  the  well  through  the  frozen  zone  there 
occurs,  on  the  one  hand,  a lowering  of  the  temperature  of  the  gas 
(oil)  due  to  the  intensive  transfer  of  heat  into  the  frozen  ground, 
and  on  the  other  hand,  the  thermal  influence  of  the  well  leads  to 
the  thawing  of  the  frozen  ground  located  next  to  the  well.  The 
thawing  boundary  shifts  with  the  passage  of  time,  which  entails  a 
change  in  the  coefficient  of  heat  transfer  from  the  gas  (oil)  in 
the  well  to  the  ground. 

Thus,  there  arises  a nonstationary  thermal  regime  in  the  well 
and  in  the  stratum.  Since  the  methodology  for  calculating  the 
thermal  regime  of  a well  in  unfrozen  ground  is  elucidated  ir,  Section 
3,  we  shall  examine  only  that  part  of  the  well  which  passes  through 
the  region  of  frozen  ground,  considering  the  temperature  of  the  gas 
(oil)  at  the  entrance  to  the  frozen  stratum  to  be  known. 

In  order  to  derive  the  fundamental  equations  describing  the 
temperature  fields  in  the  well  and  in  the  ground  surrounding  it, 
let  us  make  the  following  assumptions: 

1)  all  the  phase  transformations  in  the  frozen  ground  take  place 
at  0°C ; 

2)  the  vertical  heat  flows  in  the  ground  are  small  in  comparison 
with  the  radial  flows; 
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3)  the  vertical  heat  flow  in  the  well  due  to  heat  conduction 
and  the  heat  f*low  due  to  the  dissipation  of  mechanical  energy  are 
small  in  comparison  with  the  convective  heat  flow. 


Fig.  1. 


Key : a - thawing  boundary 

b - thawed  ground 
c - permafrost 
d - borehole 


Taking  into  account  these  assumptions  let  us  write  down  the 
equation  of  heat  conduction  for  the  frozen  and  thawed  regions  in 
the  ground.  The  heat  conduction  equation  for  the  thawed  region 
has  the  form 


dr 


**  / 


a < ' r - 'S  ft,  x) 

' n\  l 00 


(2.60) 


H9 


The  boundary  conditions  are  as  follows: 


(2.61) 


1 . V «f 


(2.62) 


The  equations  for  the  temperature  in  the  frozen  region  up  to 
the  beginning  of  thawing  are  of  the  form 


Mi  _ a j J__  + V a co 

cTl  \ ' 0}  c'r:  / 0<7<7m 


(2.63) 


The  boundary  conditions  are  as  follows: 


(2.64) 


(2.65) 


The  equation  for  the  temperature  in  the  frozen  region  after 
the  beginning  of  thawing  has  the  form 


*j--fl2(_Lie*.  + £2Ly 

dt  \r  Or  dr*  / ’ 


•S  (7,  .v)  < r < oo 

t < CO 


(2.66) 


* • *•  i 


gas  (oil)  having  in  the  initial  section  x = 0 a temperature  Tr  begins 
to  fill  the  well,  moving  with  a velocity  v. 

Thus,  the  region  of  the  change  in  the  independent  variables 
x and  t is  determined  in  the  following  manner: 


w: 


The  boundary  conditions  are  as  follows: 


T!x-o=T'0(t) . 


(2.71) 


Let  us  write  problem  (2 . 6C  )-(2 . 71 ) in  dimensionless  form. 

In  order  to  do  this  let  us  introduce  the  following  dimensionless 
variables : 


Substituting,  these  variables  into  the  equations  and  boundary 
conditions  (2. 60)-(2. 71) , we  obtain  the  following  system  of 
dimensionless  equations  determining  the  thermal  process  in  the  well 
and  the  surrounding  ground.  For  the  thawed  region  we  have: 

_ ..  ( JL  £2i  j.  1 <r<S(l,z)  (2.72) 

01  'l{  r Or  ' Or-)'  < , ^ 
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For  the  frozen  region  and  for  t — t 


sa,  /i  ae,  , tOC33  (2.75) 

at  X,-\  r Or  ^ Or*  )'  0<  t<  U 


CHAPTER  3 

THERMAL  REGIME  OF  WELLS  IN  FROZEN  GROUND 


Section  1.  Approximate  Solution  of  the  Plane  Problem  of  the 
Thawing  of  Frozen  Ground  Around  a Well 

For  high-yield  wells  the  temperature  of  gas  (oil)  varies  little 

with  time.  Indeed,  with  an  increase  in  M equation  (2.83)  approaches 

the  form 


The  solution  of  this  equation  does  not  depend  on  time,  and  for  each 
horizontal  section  of  the  stratum  it  may  be  considered  to  be  constant, 
but  varying  from  section  to  section,  and  one  may  consider  in  each 
section  the  plane  Stefan  problem.  Even  for  such  a simplified 
formulation,  in  order  to  obtain  an  exact  solution  of  this  problem, 
it  is  necessary  to  enlist  the  aid  of  numerical  methods,  which  will 
be  examined  in  the  Appendix.  Realization  of  these  methods  is  connected 
with  the  writing  of  very  complex  programs  for  a computer.  At  the 
same  time,  with  an  accuracy  sufficient  for  engineering  calculations, 
this  problem  may  be  solved  using  the  method  of  successive  approxi- 
mations described  in  Chapter  1.  Using  this  method  one  will  obtain 
a solution  directly  in  the  form  of  quadratures  for  the  case  where 
the  frozen  ground  is  at  the  thawing  temperature.  In  the  more  general 
case,  where  the  initial  temperature  of  the  frozen  ground  is  below 
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zero,  the  problem  reduces  to  a system  of  t wo  ordinary  differential 
equations,  the  solution  of  which  is  not  difficult  to  find,  in 
particular  when  one  uses  standard  programs  for  a computer. 


— . . 


Let  us  examine  in  detail  both  of  the  indicated  cases  and 
compare  the  results  with  the  exact  solution. 

Let  us  consider  the  application  of  the  indicated  method  to 
the  solution  of  problem  (2 . 72 )- (2 . 8l ) . In  order  to  do  this,  in 
conditions  (2.7*0  and  (2.77)  let  us  assume  that  T = 1.  For  the 
thawed  region  the  zeroeth  approximation  is  found  from  the  equation 


r dr  j"*"  dr- 


(3-1) 


The  boundary  conditions  are  as  follows 


eOLs  = o; 


(3.2) 


jet°> 


dfi!01!,,  - it. 


(3-3) 


Integrating  equation  (3-D  twice,  we  obtain 


«i0’  = Cx  In  r + C:. 


(3.0 
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r 


Substituting  here  the  boundary  conditions  (3-2)  and  (3-3),  we 
obtain 


rx 

a In  .S'  4-  1 


S ’ 


( ? 


1.5) 


VJe  find  the  first  approximation  from  the  equation 


, o*{l)  t <r-9<»  ^ oe;0' 

l ~7~  dr  dr- 


(3-6) 


The  boundary  conditions  are  as  follows: 


B\l)\r-S  = 0: 


(3-7) 


oe\" 

dr 


I 


= a0(,nL,. 


(3.8) 


^ (o) 

Substituting  into  (3.6)  the  expression  (3-5)  for  fcj?  , we 
obtain  the  equation 


1 \ r or  dr-  t dr 


r ^ _ g.9  (I  + ct  In  r) 


Or  ) S (a  In  S + l)'1  ' 


(3-9) 


. Oc 

Her*3  Q — — - . Integrating  the  last  eauation,  we  obtain: 


de\'} 

dr  = 2x,S  |alnS  + 


•xS 


— (r  + ar!nr  — -t)  + T-: 


(3-10) 
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9<”  = 4--I1  '"ll  4-  Cj  lnr-f  C.. 

1 *1xx6  (a  In  .S  4-  I)-  1 


(3-11) 


The  constants  and  are  obtained  from  boundary  conditions 
(3-7)  and  (3.8). 

Restricting  ourselves  to  two  approximations,  we  may  write 


% ©i0)  + e',0. 


(3.12) 


Obviously,  the  function  of  (3.12)  satisfies  boundary  conditions 
(2 . 7 3 )-( 2 . 7*0  • From  this  we  find 


<)8, 

Or 


S a 


lr.5  ~ •*  (*  In  S+l)3  Xl[2(*lnS  + ’)■ 

2a  (a  In  5 + 1)  + aa  2-~-2|r  *'-]  - — * ---  y 


(3.13) 


Analogously,  for  the  frozen  region  we  look  for  the  zeroeth 
approximation  from  the  equation 


~ — + “TP-  - 0. 


(3.1*0 


In  view  of  the  fact  that  the  condition 


(h)  — »■  as  r 

f 


oC 


which  follows  immediately  from  (2.79),  cannot  be  satisfied  by  any 


K: 
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A 


solution  of  equation  (3.14),  let  us  introduce  the  radius  of  thermal 
influence  P.(t)  [91,  over  which  we  take  the  following  conditions: 


(3.15) 


?r 


C. 


(3.16) 


Let  us  write  the  boundary  conditions  for  equation  (3.1*0  in  the 
following  form: 


o 


r*  ^ 


(3.17) 


(3-13) 


Integrating  (3.1*0  taking  into  account  (3-17)  and  (3-18),  we  obtain: 


2 6>  (3.19) 

~Tr  f (In  R-  In  S)  ' n r 


A(0)  _ llnr  — In  S) 
V-  = In/? -In  5 


(3.20) 


We  find  the  first  approximation  from  the  equation 


I J0',1'  3'0\n  \ d&l°> 


' T ‘ Sr  + Jr* 


dt 


(3.21) 


The  boundary  conditions  are  as  follows 


9 


(0 


- 0\ 


r - S 


(3-22) 


© 


(0 


= 0. 


(3-23) 


Substituting  into  (3.21)  expression  (3-20),  we  obtain  the 
equat ion 


<•  . a / 


*»  . <>  I 

r Jr  ^ 


Jr 


= 0« 


f,  ( i' 

R 

k 

, . s 

rr\— - 

K 

!:iS  — -j-  In  R 

(In  K — In  S)- 


(3.2M 


Integrating  this  equation  twice,  we  obtain: 


Jr  2 ( In  R — In  i')3  x. 


r Inr  — Sr  ) f ■ 


7T  + 


+ r ( In  S - -I"  In  R 

v K 


(3-25) 


*<•> 


4x,  (In  R — In  S)- 


+ Cj  In  r -f  C:. 


(3-26) 


By  satisfying  boundary  conditions  (3-22)  and  (3-23)  we  determine 


and  C^. 

Restricting  ourselves  to  two  approximations,  we  write: 


e2  ^ 0<o>  + q(*\ 


(3 


17) 


The  function  of  (3.27)  satisfies  the  boundary  conditions 


(3.28) 


(3-29) 


From  (3.27)  we  find,  taking  into  account  (3-19)  and  (3.25), 


+ 


60,  = W9M  [(fl-  + S-)  (In  R — In  S)  — R*  -j-  S'-\ 

dr  r=s~  4x2  RS  (In  A’  — In  S):1 

56„  [2  (In  R - In  S)  (lr.  S — In  R — I)  /e-/5-  — 1 1 ^ 9.,. 

4x,  (In  R — In  S)3  S (In  A!  — ini')' 


(9m  a 


(3-30) 


Hence  Stefan's  condition  (2.81),  taking  into  account  (3-13) 
and  (3.30),  acquires  the  form 


4 f X,9U  (2  (In  /?  — InS)  (in  S — In  R - 1)  -f-  - I] 

C>{  -lx,  (In  6!  — lnS)J 


M 

(a  In  S -f-  I)3 


2 (a  In  5 + 1 ):  — 2a  (a  In  S + 1)  4-  or  — 


2 - 2a  + a’- 1 
S- 


__  . \ K X,9„  ((R:  + S-)  (In  K - In  6')  - AT-  -4-  S-l  = 
* 4x,/v’S  (In  R — In  S)3 


X.9., 


X(a 


S (In  — In  S)  Slain  5 -i-l) 


'V\ 


(3-31) 


In  order  to  derive  the  second  equation  let  us  use  condition  (3.16) 


Taking  into  account  (3.27)  we  obtain 

R (2  (In  R - In  S)  (In  S - In  R + 1)  - 1+  SVR* | + 

4x,  (In/?  — In  5)*  (3-32) 

R 

The  problem  reduces  to  the  solution  of  a system  of  two  ordinary 
differential  equations  of  the  first  order  relative  to  the  unknown 
functions  S(t)  and  R(t)  of  (3-31)  and  (3-32).  The  boundary  conditions 
of  the  system  of  equations  (3 • 31 ) — ( 3 • 32 ) have  the  form: 


+ S 


R-  — S-  — (!<•  -f-  S-)  (In  R — In  S) 
RS 


for  ~t  tm  (3.33) 


S ^ 1 


for 


(3-3*0 


Here  R is  the  value  of  the  radius  of  thermal  influence  at  the 
m 

moment  thawing  begins.  Let  us  find  R^  and  t from  examining 
problem  ( 2 . 75 ) - ( 2 . 77 ) regarding  the  heating  of  frozen  ground  until 
the  ground  at  the  wall  of  the  borehole  begins  to,  thaw. 

Analogous  to  what  was  done  previously,  we  look  for  the  zeroeth 
approximation  as  a solution  of  equation  (3.1**)  with  boundary 
conditions 


cje(j0) 

~ir 


a(ei°>U,-ii; 


(3-35) 
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(3-36) 


& 


(o) 


r-R. 


The  solution  of  this  problem  has  the  form: 


00!, 01 


«a(0«  — 1) 


dr  r (a.  In  R + 1)  ’ 


% *'  % 


(3-37) 


q<0)  _ <9«  ~ h in  /•  4-  a.  In  R 


0 


a.  In  R -j-  1 


@M  s (3.3.8) 


Let  us  look  for  the  first  approximation  as  a solution  of 
equation  (3.21)  with  boundary  conditions: 


e£V*-o. 


(3.39) 


oe<0 

~o7" 


k=i 


a, 8 


(i) 


(3. AO) 


The  solution  of  this  problem  has  the  form: 


_‘X!KO  — Q*)('-ha1rlnr-afl-2}  , C,  . (3.^1) 

0/-  (a.  )nR  4-  ))J  ‘ r ’ 


nr 


yU>  ^ a,R  (1  — eu)r-(l  - g,  -i-  ;n,) 

4xsR  (a,  In  R -f-  t)-  ~Li  lnr  + c, 


( 3. U2) 


The  constants  of  integration  C,  and  are  determined  from  the 

1 <i 

boundary  conditions  (3-39)  and  (3.^0). 

Restricting  ourselves  to  two  approximations , we  assume 


e,  .5=.  - ei11. 


(3.43) 


Obviously,  the  function  of  (3-^3)  satisfies  boundary  conditions 
(2.77)  and  (3-15). 

Satisfying  condition  (3.16),  we  obtain,  taking  into  account 
(3-37)  and  (3.*1), 


■R.= 


\v.~  (a,  In  R — I )l 


2 (a.  In  R I )-  — 2a.  (a.  In  R — 1 ) --  a(  — • 


2 — 2a.  -j-  aj 
R1 


(3. MU) 


Hence,  satisfying  the  boundary  condition  R = 1 at  t = 0,  we 
obtain 


R, 


a.  - 2 


,xa  4^,  In  (/<■<' ct  ) 4x.a.-' In  (R-k  2 ) " l*ja3 


(3 . a5) 


Setting  equal  to  zero  the  expression  for  the  temperature  of 

the  wall  of  the  borehole  (H/  at  r = 1,  we  obtain  the  equation 

R ((2  — a.j)  In/?  4-  1 — a,  — ( 1 u-  a3  In  n?  — os,))  + (h)  — (£1-'  (3*46) 

•1  (a.  In  P -1-  1)-  x.R  In  R 4R,x.R  (a.  In  R - Ip  H p 

‘r  1-^.,  ~ a.  (l-eM)  _0- 


MM 
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Substituting  into  this  the  expression  for  R from  (3D4),  we  determine 
R , and  substituting  its  value  into  (3.45),  we  obtain  the  value 


of  t 


In  view  of  the  fact  that  the  temperature  of  frozen  rock  in 
the  majority  of  cases  is  close  to  the  thawing  temperature  (C°C  to 
-2°C),  it  is  of  interest  to  examine  the  solution  of  the  given 
problem  for  ©p  = 0.  In  this  case  equation  (3- 3D  acquires  the 


form 


^ {~  [2  ,n  'S  + ~ (a  5 + + 


, 2 — 2a  — a- 

+ «' 3^ 


1 1 _ 

J S (a  in  4-  1) 


(3.*7) 


Its  solution,  satisfying  the  boundary  condition  S = 1 at  t = 0, 


has  the  form 


. 5-  r,  ] 1 2 -2a --a- 

lxt  [ !n  (Se1  a)  j 4>'.1a-  lit  (Sis1  4) 

, S;  [2  In  (Sc1  ")-l|  a — 2 I 1 , 1 \ 


1,1) 


(3.48) 


Let  us  calculate  the  rate  of  thawing  of  the  frozen  ground  around 
a gas  welx.  Let  us  take  the  characteristics  of  the  ground  from 
data  for  the  Urengoy  deposit.  Let  us  make  the  calculation  for  the 


following  data : 


C (. ! - 0.  6 * 10  w J hour  j 3 . (3  ST  * 10  wO'  j l\cu.rj 

fr  - ; <2-  = C.H  <m  ' T-  J t0  * £Lkc^ 


' O-'O  C J ~t c ~ cO-  A CITY'S 


The  dimensionless  parameters  introduced  above  will  have  the 
values  : 

X^C.HSj  C.  flu;  X.,’  0.2 31  j C-ait. 

The  dimensionless  heat  transfer  coefficients  take  on  the  following 
values  : 

x =-  O.S'g’j  ^ * 0.4  7 3- 

The  result  of  calculating  the  rate  of  thawing  of  frozen  ground 
with  a zero  initial  temperature  using  formula  (3-48)  is  presented 
in  Fig.  2 by  a solid  line.  The  dotted  line  shows  the  exact  solution 
of  the  problem,  obtained  by  numerical  methods  using  a computer. 

From  a comparison  of  the  curves  it  is  apparent  that  the  error  in 
the  approximation  method  in  question  under  the  conditions  of  the 
given  problem  is  very  small.  For  example,  for  t = ^0  the  value 
of  S for  the  exact  solution  is  2.502,  and  for  the  approximate  one, 
2.487.  The  error  estimated  using  the  formula 

is  about  If. 

This  figure  shows,  by  means  of  a stroke-dotted  line,  the 
solution  obtained  by  the  method  of  successive  replacement  of  stationary 
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states  [5*0-  The  error  in  this  method  for  t = 40  is  about 


h cr 


Fig.  2. 

exact  solution;  2 - quasistationary  solution; 
solution  using  formula  (308). 


Let  us  return  to  the  general  case  of  a negative  initial 
temperature  of  the  frozen  ground . It  is  obvious  that  system  ( 3 * 31)  — 
(3-32)  does  not  admit  a solution  in  analytic  form.  Therefore, 
the  solution  was  obtained  on  a "Nairi"  computer  with  the  use  of 
a standard  program  for  integrating  a system  of  ordinary  differential 
equations  using  the  P.unge-Kutta  method.  In  order  to  accomplish  this, 
system  of  equations  ( 3 • 31 )- ( 3 • 32  ) was  solved  for  R and  S: 


5 = 


2ln4  + 5-l 

o A ' 


I i I R R-  St 

Mu|4ln»-r  + 2 -St-TfS 


4x.  In 


S (a  In  6'  J-  1) 

-v  k*  I r.?*  ff  / 

~JF [— - 1 +2m— I I - 


= Q 


:(n 


A., a ( 2 (a  In  5 r-  l)3  — 2a  (a  In  5 + 1)4 - a2  — - 


5 1 

2 — 2 a -u  a* 

S3 


•lx,  (a  In  5 + I )3 


(3.^9) 


i 

1= 
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R \ 


R = 


fn-g-j 


RSS 


X 


X 


(s2  - Rt  - (/?=  rS=)  In  -3- j 


, , a 

*2 In- — 


RS~- 


XjA;  (a  In  S — lj 


(aTnS-  I) 

,l2(cc  InS  -r  I): 


„ / i c * \ > » 2 — 2a  + a2  \ 
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(3-50) 


The  calculation  was  made  using  the  values  of  the  dimensionless 
parameters  indicated  above,  and  for  the  following  values  of  the 
dimensionless  initial  temperatures: 

Q - - 0.  ) ■ - 0.  ~ a JT. 


Here  Rm  had  the  values  1.432,  2.285,  and  4.220  respectively, 
and  t had  the  values  0.285,  2.0,  and  10. 31. 

The  results  of  the  calculation  are  presented  in  Fig.  3 (solid 
lines).  The  dotted  line  shows  the  results  of  the  exact  solution  of 
problem  (2.72)  — (2.  Si)  for  (§)p;  = -0.1,  obtained  by  reduction  to 
a system  of  integral  equations.  As  is  apparent  from  comparison  of 
these  results,  the  approximation  method  used  gives  very  good 
approximations  to  the  exact  solution.  Thus,  for  t = 40  the  value 
for  the  exact  solution  is  2.167,  and  for  the  approximate  solution, 
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the  value  is  2.175.  The  error  is  about  1%. 

In  order  to  obtain  an  overall  appraisal  of  the  accuracy  of  the 
method  of  successive  approximations  in  solving  Stefan's  problem, 
we  made  a comparison  of  the  solution,  found  using  the  indicated 
method,  of  the  problem  of  the  advancement  of  a plane  thawing 
boundary  (one-dimensional  case)  with  the  existing  exact  solution 
of  this  problem  [86].  The  calculation  was  made  for  the  following 
conditions: 

A ( - 0-  !(o  S j ~ 0,  1 9 b j * 0.  2 3 /j  " O.  oi'l  [ j 

" 0.  ( '}  OC  ' o<9  . 

During  the  calculations  we  restricted  ourselves  to  two  approximations. 
The  results  of  the  calculation  are  as  follows:  for  the  approxi- 

mation method  the  dimensionless  coordinate  of  the  thawing  boundary 
is  S ' ; Tor  the  exact  solution,  S = 0. 4<7Q  ft  . Thus, 

the  accuracy  of  the  method  used  here  proves  to  be  high  (about  1% 
error)  in  the  axisymmetric  as  well  as  in  the  plane  one-dimensional 
cases . 

In  conclusion  to  this  paragraph,  we  make  the  following  remark. 

In  work  [61]  it  is  confirmed  that  the  radius  of  thermal  influence 
is  proportional  to  the  radius  of  thawing  (freezing).  Fig.  t shows 
the  dependences  of  the  ratio  R/S  on  the  dimensionless  time  t for 
various  values  of  ; these  dependences  were  constructed  using 

the  results  of  the  calculation  on  a computer  of  system  ( 3 • ^9 ) - ( 3 • 50  ) . 
As  is  apparent  from  this  graph,  for  sufficiently  large  values  of 
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time,  t > 80  to  100  (i.e.  160  to  200  hours  for  our  time  scale), 
the  ratio  R/S  indeed  approaches  seme  constant  value  which  depends 
on  0^,  . Consequently,  the  confirmation  made  in  work  [61] 

is  valid  with  the  stipulation  that  the  proportionality  between  F 
and  S takes  place  for  sufficiently  large  values  of  time  and  that 
the  coefficient  of  proportionality  depends  on  ey- 


Fig.  3. 

1 - <sy=  -0.1;  2 - = -0.25;  3 - = -0.5; 

^ = -0.1  (exact  solution). 

The  idea  of  the  method  of  successive  approximations  proves  to 
be  very  fruitful  when  applied  to  the  plane  problem  of  the  thawing 
of  the  frozen  ground  around  a borehole.  As  is  apparent  from  the 
comparisons  with  the  exact  solution,  the  zeroeth  approximation 
(quasistationary  solution)  gives  an  accuracy  of  about  ^ , and  the 
first  approximation  gives  an  accuracy  on  the  order  of  1% . The 
necessity  of  the  first  approximation  is  due  not  so  much  to  this 
increase  in  accuracy,  which  may  not  be  very  significant  for  engineerin 


calculations,  as  to  the  fact  that  it  gives  the  law  governing  the 
change  in  the  radius  of  thermal  influence,  without  which  the  system 
of  equations  describing  the  motion  of  the  thawing  boundary  remains 
open.  The  solution  of  the  problem  in  question  is  obtained  either 
in  analytic  form  for  the  case  of  a zero  temperature  of  the  frozen 
ground,  or  in  the  form  of  a system  of  two  ordinary  differential 
equations  for  the  case  of  a nonzero  initial  temperature. 


Fig.  4. 

The  solution  of  the  latter  on  a computer  with  a built-in  program 
for  integrating  ordinary  differential  equations  does  not  present 
a problem.  The  fact  that  many  design  and  scientific  organizations 
have  such  computers  allows  us  to  recommend  the  use  of  the  equations 
, derived  here  for  calculating  the  rate  of  propagation  of  the  thawing 
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front  around  a borehole. 


Section  2.  Thermal  Regime  of  an  Operating  Well 

Let  us  consider  the  problem  of  determining  the  temperature 
distribution  in  a borehole  and  the  configuration  of  the  thawing 
front  at  various  moments  in  time.  Let  us  carry  out  the  solution 
separately  for  the  cases  of  oil  and  gas  wells.  In  contrast  to 
the  problem  examined  in  the  previous  paragraph,  here  the  temperature 
of  the  gas  (oil)  in  each  section  of  the  borehole  is  a function 
of  time  and  of  the  vertical  coordinate  of  the  section.  In  this 
case  the  method  developed  in  the  previous  chapter  is  not  applicable, 
since  it  assumes  that  the  temperature  of  the  gas  (oil)  in  each 


section  is  constant. 

First,  let  us  solve  the  problem  for  frozen  ground  that  is  at 
the  thawing  temperature  or  near  it.  In  addition,  we  shall  assume 
that  the  temperature  of  the  flow  at  the  entrance  to  the  stratum 
of  frozen  ground  is  constant.  In  what  follows  the  obtained  solutions 
will  be  used  for  the  more  general  case  where  one  cannot  neglect  the 
heat  flow  from  the  thawing  boundary  into  the  frozen  region. 

Let  us  rewrite  the  energy  equation  (2.63)  for  gas  (oil)  in 
the  well  in  the  following  form: 


(3-51) 


Here  the  indices  next  to  the  derivatives  indicate  which  variable  is 


held  constant  during  differentiation. 

The  rate  of  movement  of  the  thawing  boundary  is  determined 
from  Stefan's  condition: 


(3-52) 


The  temperature  field  of  the  thawed  ground  is  determined, 
assuming  a small  rate  of  movement  of  the  thawing  boundary,  from 
the  quasistationary  heat  problem: 


1 ^ _ n 

r ' Or  ' Or- 


(3-53) 


— 0‘. 


o± 1 

Or 


k=i 


= a (©i|r=i  — T). 


(3-54) 


From  this  we  find: 


ae, 

Or 


lr=I 


a T 

a In  5 - 1 ’ 


(3-55) 


je, 

or 


a T 

S (a  InS  -t-  1)  ’ 


(3-56) 
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Substituting  (3-55)  into  (3-51)  and  (3-56)  into  (3-52),  we  obtain 
the  following  system  of  equations: 


c.  (£)+ <=(£),- 


a in  5 -r  1 3 


(3-57) 


% * c 


1-5 


( 3 ■ 53  ) 


X^T 


5 (a  In  S -f  1) 


-(£)/• 


S --  i 

t-  c 


(3.59) 


(3.60) 


In  equation  (3-57)  1st  us  pass  from  the  variables  T,  z,  t to 
the  variables  T,  z,  S.  We  have: 


(§),-(f)s +(§),(§),■ 


- 


, « * J 


j 


Rough  calculations  show  that  under  the  conditions  of  the 
problem  in  question 


U%.)  Id±) 

U ft), 


« i- 


Taking  this  into  account,  and  also  using  the  expression  for 

from  (3.59),  we  obtain  equation  (3-57)  in  the  following 

form  : 


C (—)  — C (—)  — 

u 1 Us  ),  ■ (aln.S  - n ■ -(/J*  Is 


a In  S + I 


-C,; 


(3-61) 


T 


(3.62) 


Let  us  consider  the  case  of  an  oil  well  (C^  - 0).  The  equation 
of  characteristics  for  equation  (3-61)  is 


S (a  In  5 -f  I ) dS  ^ dz  _ (a  In  5 -u  1)  dT 
C\t.\%T  C2  ^7* 


(3-63) 


Hence 


SdS 


C,X 

S* 

2C,X. 


= -dT: 


= — 7~  — C. 
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For  a characteristic  passing  through  the  point  T = 1, 


we  obtain 


sl-s> 


(3-64) 


Let  us  examine  the  second  relation  of  equation  (3-63) 


S <3  In  s 4-  1)  dS  _ _££ 
Cx\xaT  C-  ' 


Using  (3.64),  we  obtain 


5 (a  In  S 4-  1)  clS 


1 


(.  S' 


s*s, 


(3.65) 


(3-66) 


The  solution  of  equations  (3.65)  and  (3-66)  gives,  in  conjunction 
vrith  (3.64),  the  desired  dependence  T(z,  s)  for  a change  in  from 
1 to  o C>  . in  the  practically  interesting  case 

KM  « L 


this  dependence  may  be  obtained  in  finite  form. 

Let  us  examine  the  function  z(S,  SQ).  From  the  boundary 


conditions  imposed  on  the  characteristic,  it  follows tha 


t ( s , Sc  1 


Let  be  such  that 


, £c)  - -1  . 


From  (3.65)  it  follows  that 


> 0 


s^s. 


From  this,  in  view  of  the  fact  that 


z & Id,  1} 


it  is  necessarily  true  that 


^ 6 R , > 


where  S,  3>  Sq.  From  an  analysis  of  expression  (3-65)  it  is 

Q)y- 

apparent  that  mon.otonically  increases  over  some  interval 

(Sp,  Sp)  from  some  finite  positive  value  to  d-  oO  . The  value  of 
S„  is  determined  from  the  relation 


Hence 


■SM  - 2 C, 


Since  due  to  the  physical  rr.eaning  cf  the  problem 
become  infinite  in  the  interval  [Sq,  S^],  S2  ^ 


2S 

V 


cannot 


Taking  into 


account  this  remark  we  obtain: 


.Sf  — S5<2C,X,; 


s — s 

O 1 ‘->0 


' 5i  -f  5( 


<c1?.1«  1. 


Consequently , 


(3-67) 


for  S 6 [SQ,  S1].  Integrating  (3.65)  taking  into  account 
(3-67),  we  obtain 


a -C. 
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S# 


K2C.it  + Si 
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In 


(C.i,  - Si  + S0  [ '2C.i,  -k Si)  <T~ 2C,j,  - Sl-S) 


- 1 + 


C.x.tv  2C,).!  + Cjj  + S) 
a In  S„  -t-  1 


(3-68) 


In 


*6-^  I 

2C,/..  |p 


Thus,  system  of  equations  ( 3 - 6U ) and  (3-68)  determines  the 
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characteristic  passing  through  the  point  T = 1,  z = 0,  S = Sq. 
Varying  the  parameter  Sp  in  these  equations  from  1 to  Ck°  , we 
obtain  the  desired  dependence  T(z,  S),  which  satisfies  equation 
( 3 - 6 1 ) with  boundary  condition  (3-62). 

Let  us  examine  the  case  of  a gas  well.  For  gas  the  effect 
of  the  inertia  term  in  the  heat  flow  equation  may  be  neglected 
[1*1,  22].  Integrating  equations  (3-61)  and  (3.62)  taking  this 
remark  into  account,  we  obtain 

f = — 1|«  C,  let  In  S — 1)|  exp  f ix  In  5 - 1)  J — 

. ' (3.69) 

- C3  (a  In  S -i-  I)J. 

Substituting  into  (3-59)  the  function  T(z,  S)  obtained  from  system 
( 3 • 6*1 )- ( 3 • 68  ) (for  oil)  or  expression  (3*69)  (for  gas),  we  find 
the  function  t(z,  S),  which  then  completes  the  solution  of  the 
problem.  In  particular,  for  gas  we  obtain 


j \a  |-  C,  (a  In  5 


S (n  In  S -I-  1 ) tlS 


"it  - n,  (a  In's  'i~)'j--C’  » 


( u ; 


And  so,  from  relations  (3.69)  and  (3.70)  one  finds  the  distribu- 
tion of  the  temperature  of  the  gas  along  the  borehole  as  a function 
of  time  T(z,  t)  and  the  configuration  of  the  thawing  boundary  at 
various  moments  in  time  S(z,  t).  The  functions  are  shown  in  Figures 
5 and  6 by  solid  lines.  The  calculation  was  made  for  the  following 
data:  M = 2.3*1  x 10^  kg/hour;  a = 0.1*lm;  = 0.03m;  clX' ~ 0.009m; 

d , = 120atn;  /V  = 0.181  keal/m-hour-deg;  c = 0.55  kcal/kg-aeg; 

■ X <X-  o 
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L = M25m.  The  remaining  data  were  taken  to  be  the  sane  as  in  the 
numerical  example  of  the  previous  paragraph.  The 

heat  transfer  coefficient,  estimated  from  criterial  dependences,  turned 
out  to  be  equal  to  p<T  = BTOkcal/m^ -hour-deg . Under  these  conditions 
the  dimensionless  parameters  acquire  the  following  values:  = 3-1^; 

C3  = 0.310;  = 0.55;  = 0.168. 

From  an  analysis  of  the  graph  of  Fig.  5 it  follows  that  the 
drop  in  the  temperature  of  the  gas  along  the*  borehole  proceeds 
most  intensively  at  the  initial  moment  in  time,  when  the  thawing 
boundary  is  right  against  the  wall  of  the  borehole.  As  the  thawing 
boundary  withdraws,  the  heat  transfer  into  the  ground,  and  con- 
sequently the  rate  of  the  drop  in  the  temperature  of  the  gas  along 
the  borehole,  decrease.  Consequently,  from  the  point  of  view  of 
the  possibility  of  hydrate  formation  in  wells  in  frozen  ground, 
the  initial  period  of  exploitation  of  the  well  is  the  most  hazardous. 

The  graph  cf  Fig.  6 illustrates  the  obvious  fact  that  the 
rate  of  thawing  of  frozen  ground  around  the  borehole  decreases 
as  daytime  is  approached,  in  connection  with  the  drop  in  the 
temperature  of  the  gas. 

In  order  to  appraise  the  accuracy  of  the  presented  method,  a 
comparison  was  carried  out  of  the  results  obtained  using  this  method 
with  the  exact  solution  of  problem  (2.72)  — (2-8^ ),  obtained  with  the 
aid  of  numerical  methods  on  a computer.  The  exact  solution  is  shown 
in  Figures  5 and  6 by  dotted  lines.  Comparison  of  the  approximate 
and  exact  solutions  shows  that  the  proposed  method  is  highly  accurate. 
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For  example,  for  t = 10  the  dimensionless  mouth  temperature  of  the 
gas  for  the  exact  solution  is  0.7765,  and  for  the  approximate  solution 
is  0.7735,  i.e.  the  error  is  about  For  large  values  of  t 

the  error  is  even  smaller.  The  error  in  determining  the  coordinate 
of  the  thawing  front  is  about  5%. 
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\ 

Fig.  5. 

approximate  solution; 
exact  solution. 


= 0;  2 - t = 10;  3 - t = 30;  H - 


Let  us  apply  the  derived  solution  to  the  case  of  a nonuniform 
frozen  stratum  with  a zero  initial  temperature.  Let  the  stratum 
consist  of  a series  of  horizontal  layers  with  different  thermophy sical 
properties,  which  are  constant  in  each  layer.  The  dimensionless 
parameters  for  the  ith  layer  we  shall  denote  by  C^,  , CxC^  , and 

/\^L^  . Let  us  denote  the  coordinate  of  the  interface  between  the 

/'j_ 


Here  T^(t)  is  the  temperature  of  the  gas  at  the  entrance  to  the 
ith  layer.  The  solution  of  problem  (3-71)  and  (3-72)  is  as 


follows : 


= t :/!  rC'3''(:‘-  lnS  + 1)|cxP(-cf‘|rn-^XT) ) - 


Cutei  lnS  -f  1);. 


(3.73) 


The  function  S(z,  t)  is  determined  from  equation  (3-59),  which  for 
the  ith  layer  is  written  in  the  form: 


i)l  J2  ' 5 (a  In  S — i) 


(3.74) 


- i 


-t*  O 


(3-75) 


Substituting  into  (3-74)  expression  (3.73),  we  obtain  an 
ordinary  differential  equation  with  the.  right  hand  side  expressed 
explicity  in  S,  s,  and  t.  Integrating  this  equation  numerically 
for  each  fixed  z,  we  obtain  the  function  S(z,  t)  in  the  ith  layer. 
Substituting  this  function  into  (3.73),  we  obtain  T(z,  t).  Recursion 
relations  for  the  temperature  of  the  gas  at  the  entrance  to  each 
layer  T^(t)  are  easily  obtained  from  (3-73),  assuming  z = Jl ^ : 


1 «"  7[.  i — j|%  T„  (I)  — C,  (7.  In -S’  [j| 

c ' p ( ~ :n ) - c ><  1 M s'  1 > ! • 


(3.76) 


Fig.  7. 

Key : a - (i  + l)th  layer 

b - borehole 
c - ith  layer 
d - 2nd  layer 
e - 1st  layer 

f - lower  boundary  of  permafrost 

Thus,  solving  successively  problem  (3-71)— (3-7^)  for  each 
layer,  one  may  determine  the  distribution  of  the  temperature  of 
the  gas  and  the  configuration  of  the  thawing  front  at  various 
moments  in  time  for  the  entire  nonuniform  stratum. 

Let  us  now  examine  the  case  where  the  initial  temperature  of 
the  frozen  ground  is  different  from  zero.  Taking  into  account 
expressions  (3.56)  and  (3-69),  we  obtain  for  a gas  well 
')e, rxT _ J_  f|  rjL  • r " 

■)r  r _ s-  i !:i  >>  1 ) S (J  a In  .S'  I 'J 

cxp(-  v.tJr*  j 

In  order  to  determine  the  heat  flow  in  the  frozen  region  let 
us  use  expression  (3-30);  then  Stefan’s  condition  (2.81)  acquires 


Adjoining  to  this  equation  (3-32),  we  obtain  a system  of 
two  ordinary  differential  equations  of  the  first  order  with  boundary 
conditions  (3*33)  and  (3-3*0.  The  variable  z appears  here  as  a 
parameter.  Thus,  this  system  may  be  integrated  for  each  value  of  z, 
i.e.  one  may  determine  the  function  S(z,  t).  Substituting  this 
relation  into  expression  (3-69),  we  obtain  T(z,  t).  System  (3-77) 
and  (3-32)  has  been  integrated  on  a "Nairi"  computer  using  a standard 
program  for  integration  using  the  Runge-Kutta  method.  In  order  to 
do  this,  equations  (3-77)  and  (3-32)  were  solved  for  S and  R.  The 
following  was  finally  obtained: 
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(3.79) 


The  boundary  conditions  are  as  follows: 


(3.80) 


System  ( 3 . 78  )- ( 3 • 79 ) was  integrated  for  the  following  values  of 
the  dimensionless  parameters:  /\f  = 0.168,  = 0.196,  = 0.271, 

= 0.55,  c 2 = 3.14,  Of  = -0.10,  = 0J72,  C3  = 0.3 1. 

The  results  of  the  calculations  are  shown  in  Figures  8 and  9, 
in  which  the  solid  lines  show  the  graphs  of  the  functions  S(z,  t)  and 
T(z,  t),  and  the  dotted  lines  show  the  results  of  the  exact  solution 
of  problem  ( 2 . 72 ) - ( 2 . 8 4 ) for  the  same  values  of  the  dimensionless 
parameters,  obtained  by  reduction  to  a system  of  integral  equations. 

As  in  the  case  of  a zero  initial  temperature,  one  may  note  the  good 
coincidence  of  the  exact  and  approximate  solutions.  The  maximum  error 
in  determining  the  temperature  of  the  gas  does  not  exceed  1*,  and 
the  maximum  error  in  determining  the  coordinate  of  the  thawing  front 
does  not  exceed  5% . 


In  conclusion  we  note  that  the  described  methodology  is  also 
applicable  for  thermal  calculations  for  boreholes  crossing  frozen 
as  well  as  unfrozen  strata.  Let  an  unfrozen  stratum  of  thickness 
be  situated  between  a frozen  area  and  a producing  bed. 


Fig.  8. 

approximate  solution; 

exact  solution. 

1 - t = 0;  2 - t = 10;  3 - t = 30;  4 - t = 50 


The  thermal  regime  of  a gas  well  in  a section  of  an  unfrozen 
stratum  is  given  by  the  formula  of  E.  B.  Chekalyuk  [19].  If  one 
calls  the  temperature  of  the  producing  bed  Tbgd,  then  the  law 
governing  the  change  in  the  temperature  at  the  entrance  to  the 
stratum  of  permafrost  will  have  the  form  [next  page]: 


*Ufr- 
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7\,v  (t)  - H.,  - C.,  In  1 1 I y.t ) 'll  H:1  - C3  In  1 1 v./l  -f"  * T 

• / , S*  e*- 
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^ V 
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Cr,  A,  are  respectively  the  specific  heat  and 
thermal  conductivity  of  the  unfrozen  ground,  and 
temperature  of  the  unfrozen  ground. 


coefficient  o 
is  the 


Fig.  9. 

approximate  solution;  exact  solution; 

1 - t = 10;  2 - t = 30;  3 - t = 5C. 
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k i.**..:.**. 


Integrating  for  the  case  of  a gas  well  equation  (3.61)  with 
the  boundary  condition 


we  obtain 


T 


^ T 6fc), 


7Y 


■!  | «~PV  (/)  . L\  (a  In  5 4-  1)1 


8<  tr^no.  e- 

(3-81) 


All  subsequent  calculations  are  carried  out  just  as  was  described 
above,  replacing  (3.69)  by  (3.81). 

Section  3.  Thermal  Regime  of  Wells  cored  in  Frosen  Rock 

In  the  presence  of  strata  of  permafrost,  the  boring  and  flushing 
of;  boreholes  using  a flushing  agent  having  a positive  temperature 
are  accompanied  by  the  thawing  of  the  frozen  ground  around  the 
borehole,  which  may  lead  to  accidental  damage  to  the  cohesiveness 
of  the  thawed  ground  [84].  On  the  other  hand, the  low  temperature  . 
of  the  strata  of  frozen  ground  causes  a significant  lowering  of  the 
temperature  of  the  flushing  fluid  and  the  formation  of  ice  plugs 
in  the  drill  pipes.  In  view  of  this,  of  much  interest  is  the 
problem  of  determining  the  temperature  distribution  in  the  flushing 
fluid  and  the  configuration  of  the  front  of  thawing  of  the  permafrost. 

In  the  present  paragraph  we  propose  an  approximation  method  for 
solving  this  problem  for  the  case  of  the  flushing  of  a borehole  bored 
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in  frozen  ground.  The  corresponding  problem  for  ordinary  (unfrozen) 
ground  was  solved  by  I.  A.  Charnyy  C 9 J - Here  his  solution,  using 
the  method  developed  in  Section  2,  is  applied  to  frozen  ground. 

As  in  Section  2,  let  us  first  give  the  solution  for  frozen  ground 
with  an  initial  zero  temperature,  and  then  let  us  apply  this  solution 
to  the  case  of  an  arbitrary  initial  temperature. 


i 


Fig.  10. 

Key : a - thawed  ground 

b - thawing  boundary 
c - frozen  ground 

Let  us  examine  the  process  of  the  flushing  of  a well  (Fig.  10) 


bored  in  uniform  frozen  ground  with  an  initial  temperature  of  0°C. 

After  C 9 J , the  heat  flow  equation  for  the  flushing  fluid  in  the 
central  pipe  and  in  the  annular  space  is  written  as  follows  [next  page]: 


The  heat  conduction  equation  for  the  thawed  ground,  in  accordance 
with  the  assumptions  made  in  Section  4 of  Chapter  i,  is  written 
as  follows: 


(3-85) 


The  boundary  conditions  are  as  fellows: 

«.k  - ” 0: 


(3.86) 


(\J  3 


The  rate  of  advancement  of  the  ground  thawing  boundary  is  determined 
from  Stefan's  condition: 


i)S 


Or 


(3. 38) 


Let  us  introduce  the  dimensionless  quantities 


T — • r = _!_•  (.)  h,  . s_,  ^ 

' 7'„x‘  / ’ >'  ’ 1 7'„v ' « a, 


fc;  • Al 


~r 


e ra,  r i>  c. 


In  these  variables  equations  (3-82)— (3- 88)  acquire  the  form: 


21*1  - - iiiiu  bl-  ,,  ,T  _ r )■ 

, , '-I  I'  1 ' 


(3.89) 


I ~l)  « — •»,  (7*,  - r.); 

* I , 


(3.90) 


i;i"i  _ / I 1)1-1, 

1,1  y'[  [ r ilr  ’ ilr - 


(3.91) 


*- 

c 

t 

f 

■ f> 
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In  view  of  the  fact  that  the  rate  of  advancement  of  the  ground 

thawing  boundary  is  small,  one  may  consider  that  the  distribution 
* 

of  the  temperatures  in  the  thawed  region  differs  only  slightly  from 
the  stationary  distribution.  Taking  into  account  this  remark 
we  obtain 

I 


<£ij af;  oH,  I ar, 

- |r»s  ' SlainS  — i)  ’ Or  |,  = , a In  S ~ ] 4 

I 1 

Substituting  relations  (3*9*0  into  (3-89),  (3-90),  and  (3-92), 


we  obtain: 


*9*0 

1 


Mr. 

•S'  I 'J.  In  .s  ■ I)' 


(3.97) 


Thus,  the  problem  reduces  to  the  solution  of  a system  of  the  three 
differential  equations  ( 3 . 95  )- ( 3 . 97  ) with  boundary  and  initial 
conditions  (3-93)  - 

In  system  ( 3 . 95 )- ( 3 • 96  ) let  us  pass  from  the  variables 
(Tlt  T2,  z,  t)  to  the  variables  (T.^  T2,  z,  S),  taking  into  account 
that 


Rough  calculations  show  that  for  boreholes  having  a sufficiently 
large  diameter  ( 190mm),  -fc-  is  sma • Taking 

this  into  account,  equations  (3-95)  and  (3.96)  may  be  rewritten  in 
the  form: 


*,t2 

a In  5 1 


■Jl-T,)-. 


(3.98) 


’ — «i  (T,  - T.,). 


(3-99) 
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Eliminating  from  this  system  , T^,  we  obtain 


L (* 

f t 

; . 

I 

* 5 , 


— a-a' r 1 0 

\ <)z-  Js  a la  S ;•  t ( oz  j5  otln-S-i-l 


(3-100) 


The  boundary  conditions  are  as  follows: 


- 1- 


^ --o 


(3.101) 


9~r, 


9 Hr 


- 0. 


(3.102) 

The  latter  condition  is  obtained  at  once  from  equation  (3*99)  and 
boundary  conditions  (3.93). 

Solving  for  each  fixed  S the  boundary-value  problem  (-3-100)  — 
(3.102),  we  obtain 


where 


r >e  * — /,/• 


(3.103) 


1 ' 

| 1 (<x  In  .V  -|-  1 )'  1 « 1 11  -S’ 4 l 1 


“ 3 (a  In  .S  l) 
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Analogously,  for  T,  we  obtain  the  equation 


( Oi1  J s i In  -S'  r 1 V l>!  )s  '*!»•!>  I 


- 0. 


(3.104) 


The  boundary  conditions  are  as  follows: 


r.Ui  - 


(3.105) 


1ZL* 

* 2-1 


q;r, 

a In  6'  -r  I 


(3.106) 


The  solution  of  problem  ( 3 . 104 )- ( 3 . 106  ) for  each  fixed  S has 
the  form 


T, 


rxcr'  - /•,«'> 


(3.107) 


Thus,  we  have  obtained  T2(z,  S).  Sibstituting  (3-107)  into 
(3.97),  we  obtain  an  ordinary  differential  equation,  the  right 
hand  side  of  which  is  a known  function  of  z and  S.  Solving  this  equa- 
tion numerically  for  each  fixed  z,  we  obtain  the  function  S(z,  t). 
Substituting  this  expression  into  (3-93)  and  (3.107),  we  obtain 
T^(z,  t)  and  T2(z,  t). 

Let  us  examine  the  case  of  the  flushing  of  a borehole  that  is 
bored  in  frozen  ground  with  an  initial  temperature  different  from 


zero.  Taking  into  account  expressions  (3.9^)  and  (3.107)  we  obtain 
for  the  thawed  region 


(>-*i 

Or 


*T; 

J>  (i  In  S I ) 


i (a  I n 6' t-  1 ) 


1 — r,e 


(3-108) 


In  order  to  find  the  heat  flow  into  the  frozen  region  let’ us  use 
expression  (3.30).  Then  Stefan's  condition  (2.8l)  acquires  the 
form 


•,  j ;';e"  [2  ln  s ln  Rc  + s-  ~ 


lx.  In1  -$r 


- 1 


+ R 


.*>•)  In  ~§~  — R: 


I y...RS  In1  — 


a-9>, 
5 In  ~ 


NlllnS  I) 


&lt  « <3)c 


KK 


(3.109) 


r .<  ■ - r.r 
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Adding  here  equation  (3.32),  we  obtain  a system  of  two  ordinary 
differential  equations  relative  to  the  functions  S(z,  t)  and  R(z,  t). 
As  initial  conditions  we  may  use  conditions  (3-33)  and  (3.3^). 
Solving  this  system  for  each  fixed  z,  we  obtain  the  desired  function 
S(z,  t),  and,  substituting  this  function  into  expressions  (3.103) 
and  (3.107),  we  obtain  T^z,  t)  and  T2(z,  t),  which  then  completes 
the  solution  of  the  problem. 

Let  us  calculate  the  temperature  of  the  flushing  fluid  during 
the  flushing  of  a well  in  frozen  ground  with  a zero  initial  tem- 
perature. The  calculation  will  be  made  for  the  following  beginning 


data:  Tentrance  = 20°C;  L = 440m;  a = 0.095m;  b = 0.052m;  7^  = 

= 0.8  x 10”^m2/sec;  ^ j - 0.147  x 10'^kcal/m-sec-deg;  k1_2  = 0.396 

kcal/m2-sec-deg;  pC  ~ 0.760  kcal/m2-sec-deg;  ^ = 0.428  x 10  ^ 

kcal/m-deg-sec ; c = 1 kcal/dg-deg;  tg  = 2 hours;  G = l8kg/sec.. 


z 


Fig.  11. 


The  results  of  the  calculation  are  presented  in  Figures  11 
and  12.  The  solid  lines  show  the  distributions  of  the  temperatures 
T1  and  T^,  and  also  the  configuration  of  the  thawing  front  S,  calcu- 
lated for  various  moments  in  time.  The  dotted  lines  show  the  results 
of  the  exact  solution  of  problem  ( 3 . 89 M3 . 93 ) , obtained  using 
numerical  methods.  Comparison  of  the  approximate  and  exact  solutions 
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shows  that  the  applied  method  is  highly  accurate.  Thus,  for  t = A 
the  temperature  of  the  flushing  fluid  at  the  exit  from  the  annular 
space  for  the  exact  solution  is  0.919,  and  for  the  approximate 
solution  the  value  is  0.931-  Consequently,  the  error  in  determining 
the  coordinate  of  the  thawing  front  is  about  1.2%. 

Section  4.  Calculation  of  the  Rate  of  Freezing  of  Previously 
Thawed  Ground  After  Shutdov/n  of  a Well 

After  cessation  of  exploitation  or  flushing  of  a well,  the 
liquid  remaining  in  the  well  rapidly  acquires  the  temperature  of 
the  surrounding  ground  and  the  previously  thawed  ground  begins  to 
freeze.  The  result  of  this  process  may  be  the  appearance  of 
significant  normal  stresses  on  the  wall  of  the  borehole. 

Frozen  ground  is  a non-Newtonian,  elastoviscous  medium,  i.e. 
the  stresses  in  it  depend  on  the  rate  of  deformation  and  on  time. 
Since  the  deformations  in  the  given  case  depend  on  the  volume  of 
the  freezing  ground,  the  rate  of  freezing  has  a considerable  effect 
on  the  magnitude  of  the  stresses  that  arise  in  the  ground  and 
on  the  wall  of  the  borehole.  The  problem  of  the  rate  of  freezing 
of  the  ground  around  a borehole  thus  acquires  great  practical 
importance . 


The  temperature  field  in  the  thawed  ground  is  descr 
equation  (2.72)  with  boundary  conditions  (2.73)  and  - 


Lbed  by 

- 0 

r-  L 

i.e.  we  assume  that  the  temperature  of  the  liquid  rapidly  drops  down 


to  the  temperature  of  the  surrounding  ground,  in  connection  with 
which  the  heat  flow  at  the  wall  of  the  borehole  will  be  equal  to  zero. 
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Fig.  12. 

The  temperature  field  in  the  frozen  ground  is  described  by- 
equation  (2.78)  with  boundary  condition  (2.80)  and  initial  condition 
(2.79),  where  by  tm  will  be  understood  the  time  at  which  freezing 
begins,  and  by  will  be  understood  the  temperature  field  in 

the  frozen  ground  at  this  time.  Here  it  is  necessary  to  add  a condition 
at  the  freezing  boundary  [(2. 81)]. 

As  rough  calculations  show,  the  temperature  of  the  thawed  , 

ground  after  a very  short  time  interval  (in  comparison  with  the 
total  freezing  time)  becomes  close  to  0°C,  i.e.  during  practically 
all  of  the  freezing  time  the  flow  from  the  thawed  region  is  small 
and  may  be  neglected  when  compared  with  the  heat  flow  in  the  direction 
of  the  frozen  ground.  An  analogous  assumption  is  made  in  work  [5^]. 
Thus,  the  condition  at  the  freezing  boundary  may  be  written  as 


¥ 


Solving  the  heat  problem  in  the  frozen  ground  using  Ifr'  method 
presented  in  Section  1,  we  obtain,  taking  into  account  (3-110), 

_ *«.[<«■  + *>'•■' f -*■'■*]  «ffl  « & 

dr  4x.RS  in3  -j- 

f . R i.  S \ , .1  (3-111) 


Adding  to  this  equation  (3-32),  we  obtain  a system  of  two  ordinary 

differential  equations  of  the  first  order  relative  to  the  functions 

• • 

S(t)  and  R ( t ) . Solving  this  system  for  S and  R,  we  finally  obtain: 

/ R „ R:  S'-  \ 

. ,a  / » « \ XtH«(4ln*T+2“  S*  -«*) 

S,^(  ..n-f  + Jr-.)  ^ 

« i,  , it  u S 


^-0 1 v 

>.,9u(/?3-S3-2S*lii-f-]  4*,  In1  -j  1 | 

#■ E* +_ Tt H 

|\  4X«ln~  ‘ 


(3-112) 


(3-113) 


R\‘-‘,n  ~ Rm’ 


(3-11M 
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Here  the  values  for  R and  S are  obtained  from  the  solution  of 

m m 

the  problem  of  thawing. 


s i 


Fig.  13. 

1 - ©F=  -0.25;  2 - fij-*  -0.1. 

System  (3  • H2)-(3  • 113  ) was  integrated  using  a standard  program 

on  a "Nairi"  computer  for  the  values  of  the  dimensionless  parameters 

cited  above.  The  results  of  the  calculation  are  shown  in  Fig.  13- 

From  these  results  one  may  conclude  that  the  freezing  of  ground 

that  had  thawed  to  some  radius  S requires  much  more  time  than  the 

m 

thawing  to  this  radius.  For  example,  for  (£$),  = -0.1  and  S^  = 3 
the  thawing  time  is  168  hours,  and  the  refreezing  time  is  900  hours. 
This  may  easily  be  explained  from  the  following  considerations. 

The  flushing  and  exploitation  of  wells  is  an  intensive  thermal 
process,  connected  with  the  effect  of  large  amounts  of  warm  fluid  on 
frozen  ground.  In  this  connection,  the  heat  flow  at  the  wall  of  the 
borehole  is  great,  and  the  rate  of  thawing  is  sufficiently  great. 


At  the  same  time,  the  process  of  refreezing  is  much  less  intensive, 
since  the  temperature  of  the  frozen  ground  usually  differs  little 
from  0°C,and  the  heat  flow  toward  the  frozen  ground,  which  determines 
the  rate  of  freezing,  is  small. 


APFENDIX  1 


1.  Reduction  of  the  Problem  of  the  Heat  Exchange  Between  a Well 
and  Frozen  Ground  to  a System  of  Integral  Equations 

In  order  to  appraise  the  accuracy  of  the  above-developed 
approximation  methods  for  calculating  the  thermal  regime  of  wells 
in  a region  of  permafrost,  it  is  necessary  to  compare  the  results 
obtained  using  these  methods  with  the  exact  solution  of  problem 
(2 . 72  )- (2 . 84  ) . In  order  to  obtain  the  latter  let  us  use  the  method 
of  reducing  the  axisymmetric  Stefan  problem  to  a system  of  integral 
equations  [21].  Work  [21]  examines  the  one-dimensional  axisymmetric 
problem  of  the  melting  of  a cylinder  of  infinite  length,  the 
temperature  of  the  surface  of  which  varies  according  to  an  arbitrary 
law.  In  cur  case  we  are  concerned  with  two-dimensional  problems 
in  a semi-infinite  region  with  a cylindrical  hole,  in  which  the 
temperature  varies  with  time  and  with  height.  However,  using  the 
methodology  developed  in  [21]  and  [86],  it  proves  to  be  possible  to 
reduce  this  problem  to  a system  of  Integral  and  ordinary  differential 
equations . 

Let  us  examine  the  region  EAEF  (Fig.  14),  bounded  by  the  char- 
acteristics AB  and  EF  and  the  curves  AE:  r = r^(t)  and  BF:  r = r2(t). 
Let  us  consider  sufficiently  smooth  functions  f such  that 


According  to  Green's  formula,  for  the  region  PABQ  we  have 


Hence 


\ \r**dr  + 

PAUQl* 

-AM** +*)-**{**+#) 


r T 


•ft 


~n J>5-  drdt-  fj'[nf8(<f)-r(f‘»(ij,)|drd/-0. 


. J r<r>\dr=  j nfi\-dr  ■-  j -f-  a-r  - «p  )<*/]■  (A.lj 

f’Q  A/3  BQ'i'PA1" 

Let  <p  (r,t)  * ©(r,  -t)  be  some  solution  of  the  heat  conduction 
equation.  Let  us  examine  the  function  of  an  instantaneous  cylindrical 
source  of  heat 

„ mp(-£(7=T)  , ( rl  \ 

En  (r.  «,  x (/  - t))  = ’* ( 2x17^7)  )• 


For  this  function  we  have 


8(£fl)=0  over  r o.vJL  ~t j 

OT(£o)=0  over  ^ 'tJ>% 

Let  us  examine  the  region  PABQ  in  the  coordinates  ^ , 3?  and 
the  points  M(r,  t)  and  M1(r,  t + h).  Let  us  assume  that 

rp  = (-)(j.  t);  i|’  = £0 (r.  |.  x (t+h— t)). 


0- 


The  last  limit  is  substantiated  in  [86].  Hence 


0 (r,  I)  = \ •(•)  (6,  /)  £„  (r,  x (t  - t))  d£  + 


+ \ [|0(s.  5. 


(A. 2) 


Let  us  examine  the  application  of  the  general  relation  (A. 2) 


to  problem  ( 2 . 72  )-( 2 . 8l ) . For  a thawed  zone  region  PABQ  acquires 


the  form  shown  in  Fig.  15-  Relation  (A. 2)  for  the  indicated 


region  is  written  as  follows  (here  and  below  the  index  z is  omitted) 


0i  (r.  t)  = \ ge,  (|,  t)  £0  (r.  $,  x,  (/  - T))  dl  + 


I3Q 


M BQ  ' ' 

•m 

" 1 |*l(£n('.  I.  X,(?-T)))^'  -/( T)2J°  ]dt- 

y *•  51  S-i  "s  |-i. 


(A. 3) 


+ 1 *iS(t)  £„(r,  S(t),  xt(f  — t))-  ^’l  dr. 


Here 


/(0«e.(i.  0- 


(A.U) 


For  the  frozen  region  of  problem  (2 .72 )-(2 . 8l ) the  points  B and 


Q of  region  PABQ  are  removed  to  infinity.  In  order  to  show  that  it 
is  possible  to  apply  relation  (A. 2)  to  a region  of  this  form,  let  us 


□ 
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n 

This  takes  place  for  any  t > tm,  since  — ^ - j _ ^ ^ as  R — 0 
(the  perturbations  of  the  temperature  field  at  infinity  die  down)  and 
the  exponential  factor  goes  to  - oo  as  R — > . 

We  note  that  this  integral  is  improper,  with  a removable 
singularity  at  the  point  'tr  = t . As  T— ^ t the  integrand  goes 
to  zero: 


t)^|  (/t=-  f x.K«h  </?.  t) 

' r:  f i?*  \ 

,Apl~  lx  <./  — t;  / / A'  i / '><  ) ' v 

> >x„  ((  — T)  l Ux.  (/  — T)  0 \ 'Jx,  (/  --  T)  / Tl 

- /.  "I’ ( “ I ■ • 

'•(sTi^To )) I ■H** **•  - 

*■11 

X [(r  --/?)( I -i-  0 f -^))]  d r — ^ 0. 


The  region  PABQ  for  the  frozen  region  may  be  obtained  from 


region  SR  for  R 


. Hence 


0,  (r,  0 = j’  ?0;  (j,  *m)  £„  (r,  s.  x,  (/  - /,„))  iiS  - 
1 
t 

- j y*S  (r)  £„  (r,  S (r),  x,  (/  - t))  p dr. 

i...  i S(T) 


(A. 6) 


Introducing  the  notation 


OH, 

«.i'i 


MO-S* 


L . * . . H • U 


we  obtain 


i 

0,  (r.  0 = fx.,5  (t)  £„  (r , S (x),  xt  (/  - x)>  v^dx  - 

f / 

■ ] x.,£a(r.  1,  x,  (t  — i))  vt(x)dT~-  [x,/(t)^|  dT 

<m  <«  * U~' 

= S,  - S.  n Sr 


(A. 7) 


0,  (r.  t,  =•  j &«0  (')  5,  X-  (/- 1 „))  Cil  - 

i 

— | x,S  (x)  E0  ( r , S (r),  xs  (<  — > ■)  v3  (t)  dr  = r . — 5V 


(A. 8) 


Let  us  verify  the  possibility  of  differentiating  the  improper 
integrals  through  with  respect  to  the  parameter  r.  Integral 
converges  at  any  point  in  the  interval  1 ^ r S(t),  and  the 

result  of  formal  differentiation  of  this  integral  under  the  integral 
sign  converges  uniformly  on  any  interval  1 4s  4.  r 4-  r^  ^ S(t). 
Integrals  S2  and  converge  at  any  point  in  '.he  interval  1 4=  r 4s.  S(t 
and  the  results  of  formal  differentiation  of  these  integrals  under 
the  integral  sign  converge  uniformly  on  any  interval  1 4z  r^  4 r 4s 
4s  r^  4s-  S(t).  Integral  converges,  and  the  result  of  formal 
differentiation  of  this  integral  under  the  integral  sign  converges 
uniformly  on  any  interval  S(t)  4 ^ ^ r 4 r?  4 o£)  . Integral 
converges,  and  the  result  of  formal  differentiation  of  this  integral 
under  the  integral  sign  converges  uniformly  on  any  interval  S(t) 
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r1  ^ r r2  ^ cO  ■ 


Thus,  on  any  interval  1 ^ r^  4=  v £ 


?2"^-  s ^ » 


t* 

f 

« 

* *• 
li  v. 


— - 1 Xji*  (T)  :•,  (T)  E0  (r,  S (t).  x,  (/  — r))  d t — 


/ ; 
- \ y.^Jx)—  f«(r>  !•  x,(/-T))dT  - 


Taking  into  account  the  fact  that 


<W:*  ^ /T  ( 

*l  r5T* 


the  last  integral  may  be  transformed  as  follows: 

/ { 

■'  f x./ <T) t/T  = \f(x)dEx  = f(l)-Ex(r,  1,  0)  — 

— / (/m)  fi  (r,  1.  x,  (/  — r,,,))  — \ C,  (r,  1,  xt  (/  — t))/(t)  dr. 

•n, 

Let  us  show  that  the  first  two  terms  are  equal  to  zero 


Ei  (r.  1,0)—  !im- 


exp 


r*  I 


A|V  <*■' 


»»  liin 


«,<  '.(rfer) 

"i’(  •Tw1)"r(-^r)f'  ' 

= lirn ' = 

rr)..  cm, 


I .0  I'.’x,/)1  •'  | .'.V 


— 0 irpn  r I/-,,  r . 


112 


On  the  other  hand,  f(tm)  = 0,  since  this  is  the  temperature 
at  the  wall  at  the  moment  thawing  begins.  We  finally  obtain 
for  r & [r1,  r?] 


1!F  = J (t)  Vi  (t)  ± e0  (r,  S (r).  z,  (t  - t))  dx  - 

‘m 

I 

- J (T)  -J?  £»  (r.  1.  Xt(l  - r))dx  — 

t 

J (r,  1,  v.j  (/  — x))  / (t)  dx . 


(A. 9) 


Differentiating  (A. 8)  with  respect  to  r on  any  interval 


S(t)  Z r^  4=  r £ T 2 > we  obtain 


5 « f |0O  (s)  4r  £„  (r.  s.  (l  - U)  dl  - 


(*■  ‘'O 


— J x.S  (t)  fj  (t)  -jp  £0  (r,  S (x),  x.  (/  — r))dt. 


Taking  into  account  that 


<)r 


££i_£1 

(/£  • 1 


t 

h* 

. ' «t 


let  us  transform  the  first  integral  in  the  following  manner: 

CD 

7 -Is©. ($)■!■£, (r.i.  x„(C-/m))d?  = 

*0  CN» 

= ,f  se#  (g)  Jjjjji  <fg  - J'  0O  (g)  Etdl  = |6n  (;)  £,|r  + 

CO 

- J £,</(g«,(g))- Je.tgjE.dg. 
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Substituting  the  values  of  the  limits  of  integration  into  the  first 
term  and  taking  into  account  that 


I'm  C,  [r,  l,  x.  (/—  ir 


.)>  - 0 


and 

&c  0-)  * 


we  obtain 


' = J (<->„ <s ) ! !<-h  (?))  d6  - [ e„  (D  E,d\ 
' i 

\ 

= \ (r,  c,  x...  (/-/,„))  i©o  (?)  dfi. 


Hence  for  r [ r7  , r2]  we  have 


- J £\(h  Si  x.(/-/m))s0o (s )dl  - 

t 

— j x.,5  (r)  i\  (t)  —■  En(r,  S (r),  x,  (f  — x))dx. 

t,„ 


(A. 10) 


Let  us  consider 


.. 

Ian  -r— 

r.  I ^ 


Integral  in  expression  (A. 7)  is  a differentiable  function  of 
the  parameter  r on  any  interval  1 h r ^ v ^ 


r2  ^ S(t ) . Hence 


Let  us  consider  integral  S^.  Using  the  integral  representation  of 
the  Eessel  function  of  an  imaginary  argument,  we  transform  this  integral 
in  the  following  manner: 


A 


From  this  we  obtain  (Fig.  17) 


Fig.  17. 

Here  ^ is  the  unit  circle.  Thus,  integral  S£  is  in  the  form  of 
the  two-dimensional  thermal  potential  of  a simple  layer  [87]. 

U.-.ing  the  theorem  regarding  the  discontinuities  of  the  normal 
derivative  of  the  thermal  potential  of  a simple  layer  [88],  we  obtain 


hm  — ? = 

, .1  Or 


I.  II)  , 

2 of 


r«l 


The  last  integral  in  expression  (A. 9)  uniformly  converges  in  r 
on  an  Interval  which  includes  the  point  r = 1.  Let  us  show  this, 
Let  us  examine  the  behavior  of  the  integrand  as  T —>•  t . 


/ (t)  £j  (r,  !,  x.({  — t))  =/( t) 


exp 


\ '1*1  U-T)  v . 

— X 

— ri 


'•(sri7-«)~,’(,) 


/ r*  I 1 \ / r \ 

(v-  TT7777— 77 ) ( ^777^7- ,) 

-*i  0 I :!.i  - — — - 

> l’Xl  (t  — T) 


(r  ]i-  \ 


( , , „ (.  + » (***))< 


: / (t) 


|'  i:\ry,  (/  — rt 


npn  r " 1 1,  r2]. 
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It  is  obvious  that  the  integral  of  the  last  function  (maj orante ) 
converges  (f(t-)  is  bounded),  and  consequently  the  integral  in 
question  converges  uniformly  on  any  interval  Cl,  r^]-  Hence 


3 i rn  j /(t) (r,  I,  x,  (/- t))rfT  j / (t)  £,  (1,  1,  t)) di. 


We  finally  obtain 


lim  = ;-„(/)  = j x,.S'  (t)  (t)  ~E0(\,  S (z),  xl  {t—r))dz  ~ 


r * I 


4.  - ]’  xltf„  (r)  ■—  £u  (I,  1,  x.j  (/  -x))  dr  - 


— f /(t)  £j  (1,  1,  x,  (/— T))dt. 


Let  us  examine 


A 


?e, 


OUr\ 


as  A 


S>C<=). 


PA 


Analogous  to  what  was  done  above,  one  may  show  that  the  integral 
in  (A. 7)  may  be  put  into  the  form  (Fig.  18) 


i f exp(  ■)*,  (/— T)  ) 

) ^(x)-  \--~-dydz. 


According  to  the  theorem  regarding  the  discontinuity  of  the  normal 
derivative  of  the  thermal  potential  of  a simple  layer  we 
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» H,  J u **  ■ .4 


obtain 


I ini 

r » >'W) 


6S  { 
7*7 


r.  U) 


+ ‘— ' 
■2  ' rtr 


r — S(l) 


It  is  easy  to  show  that  the  second  and  third  integrals  in 
expression  (A. 9)  uniformly  converge  in  r on  any  interval  1 r^  ^ 
^ r kl  s(t),  and  consequently  are  continuous  in  r at  the  point 
r = S(t).  Taking  this  into  account  we  obtain 


lim  97  - yi  (*)  = ■"  1 '-is  (t)  i'i(r) • 4-  i'0  ■ 

in 

. t 

X (5  (/),  S (r),  Xj  (t—x))dx  — j x1»0  (x)  E„  \ 

t 

' is  (t)>  I,  y-t  (t— T))  dx — [ / (x)  £,  (S  (/),  1,  (t—x))dx. 


Consider 


, • ()W„ 

lun  — 

r .S(/l 


As  was  done  above,  one  may  show  that  integral  in  expression 
(A. 8)  may  be  put  into  the  form  of  the  thermal  potential  of  a simple 
layer  (Fig.  19) 


S3 


lx,  (/ — T) 


•».T  (/  - X| 


dydx. 
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Fig.  18. 


Fig.  19. 


Hence 


, • "S-. 

lim  — 

r ..MO  "r 


(t ) 
2 


The  first  integral  in  expression  (A. 10)  uniformly  converges 
in  r on  any  interval  S(t)  r r^  <o©  , and  consequently  is 

continuous  in  r at  the  point  r = S(t).  Hence  we  finally  obtain 


lim 


TF  = (0  = ] (5  (0. 1 X,  (/-/.„))  ge;  (?)  di  H-  ^ 

I * 

/ 

— J x,S  (t)  i’2  (r)  — £0  (5  (7),  5 (r),  x{  (/ — t))  7t. 


Thus,  problem  (2 . 72  )-( 2 . 80  ) reduces  to  a system  of  integral 
equations  of  the  Volterra  type: 


v0(t,  z)  = 2 fxj  (z)  S (r,  2)  0l  (t,  2)|t£0(1.  S(/— r),  x^)  (/— r))dr_ 

t t 

-2  jy.l(z)v0(T,2)-£rE<t(l,  I,  Xt  (2)  (/  — t))  dx  — 2 f//( r,2)X 

^ Vi 

X £i(l,  1.  X[  (2)  (/ — t))  lit; 


(A. 11) 


119 


■S  (/,  z)  = I r ] i — - t (r,  z)  -r  /•;  (z)  v2  It,  r)|  dz.  ( A . 14  ) 


and  the  expressior  for  the  temperature  of  the  wall  of  the  borehole, 
according  to  (2.74): 

f(t,z)  = '^~-TU,z).  (A.  15) 


i> 

e 
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The  system  is  completed  by  the  equation  for  the  temperature  of 
the  gas,  which  may  be  obtained  from  (2.83)  , if  one  neglects,  as 
is  usually  done  for  gas,  the  inertia  term 

TU,  7)»  1 \ (A. 16) 

j 

System  of  equations  (A. 11) - (A. 16)  was  solved  using  the 
method  of  successive  approximations  on  a BESM-3M  computer.  In 
order  to  do  this,  a program  was  written  in  the  ALGOL  language  for 
the  o<-  -translator.  Let  us  clarify  the  basic  principles  on  which 
the  program  works.  In  each  horizontal  section  z = z^,  beginning 
with  the  section  z = 0,  plane  problem  (A. 11) - (A. 15)  was  solved. 
Using  the  function  v^  obtained  from  this  solution,  from  equation 
(A. 16)  the  function  T was  determined  for  the  subsequent  section 
z.  , = z.  + A*?:  . This  function  was  substituted  into  system 
(A. 11) - (A. 15) , and  the  values  of  all  the  desired  functions  for 
this  section  were  determined.  The  procedure  was  then  repeated. 

In  order  to  solve  plane  problem  (A. 11) - (A. 15) , n points  on  the  time 
axis  were  singled  out.  The  coordinate  of  the  mth  point  was 
denoted. by  D[m],  where  m = 1,  2,  ...,  n.  The  values  of  the  unknown 
functions  v^ , v^,  v^ , s,  f,  and  T at  these  points  were  denoted 
respectively  by  vOj  1 [m,  p]  , vljl  [m,  p]  , v2jl  [m,  p]  , sj  [m,  p]  , 
fj  [m,  p]  , Tgj[m,  p]  (previous  approximation)  and  v0j[m,  p]  , 
vljfm,  p]  , v2j[m,  p]  , sjl[m,  p]  , fjl[m,  p]  (next  approximation). 
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In  order  to  construct  the  function  over  the  entire  axis 
from  its  values  at  n points,  the  interpolation  procedure  w was 
introduced,  which  to  each  of  the  discrete  functions  vOjl  [m , p]  , 
vljl  (m,  p] , . ..,  determined  at  the  n points,  put  in  correspondence 
the  continuous  functions  w[y,  vOjl]  ...,  determined  in  the  following 
way.  At  the  reference  points  D[m]  the  values  of  the  functions 
w [y , ...]  coincided  with  the  values  of  the  corresponding  discrete 
functions,  and  at  intermediate  points  the  values  were  determined 
by  linear  interpolation.  For  example,  the  function  w [y,  vOjl] 
is  a continuous  function,  defined  in  the  following  manner:  for 

y Q_  [D  [m]  , D [m+1]  ] , 

*’['/.  t-0/I]  = (fO/I  [/«  + [,  p]X U/-D [.n])  + 

+ vQjt[m,  p\X(V[iii+l\—y))f{D[m+\]—D[m]\. 

Using  the  functions  w[y,  ...],  f or  each  mth  point  were  constructed 

t 

the  functions  zO , zl,  z2  , which  are  the  integrands  in  the  integrals 
over  the  time  axis  in  equations  (A. 11)  , (A. 12)  , and  (A. 13)  respectively, 
and  z3,  for  the  volume  integral  in  equation  (A. 13).  The  functions 
vOjfm,  p]  , vlj[m,  p]  were  determined  by  integrating  the  functions 
zO  and  zl  respectively  within  the  limits  D[l]  to  D[m].  The  function 
v2j  [m,  p]  was  determined  analogously,  with  this  difference:  tc  the 

integral  of  the  function  z2  within  the  limits  D[l]  to  D[m]  was 
added  the  volume  integral  of  the  function  z3. 

The  time  integrals  in  equations  (A. 11) - (A. 13)  are  improper, 
with  a singularity  at  the  upper  limit  of  the  form  1/  ft  • 
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Correspondingly,  the  functions  zO  , zl,  z2  also  have  at  the 
points  D[m]  singularities  of  this  form.  Taking  this  into  account, 
the  region  of  integration  was  broken  up  into  two  intervals:  [D[l]  , 

D [m]  - £.  1 [m]  ] and  [D[m]  - £l[m],  D[m]],  where  £ 1 [m]  is  a 

value  that  is  so  small  that  in  the  last  interval  the  integrand 
may,  with  a sufficient  degree  of  accuracy,  be  approximated  by 
a function  of  the  form 

In  the  first  interval  the  integral  was  calculated  using  the 
standard  procedure  "simps;"  the  integral  in  the  second  interval 
was  calculated  analytically  from  the  approximating  function 

, the  coefficients  in  which  were  determined  from 
the  values  of  the  integrand  at  the  points  D[m]  - £,l[m]  and  D[m]  - 

tl[n]/2. 

Thus,  given  some  first  approximation,  i.e.  a set  of  discrete 
functions  vOjlfm,  p]  , vljlfm,  p] , . ..,  successively , for  each  m, 
we  obtain  using  these  functions  a new  set  of  functions  vOj[m,  p]  , 
vlj[m,  p] , ...  . The  second  approximation  is  taken  in  the  form 

: ’ I ["'■  P\  ■!-  :P;  ! "•  n\ 

■>  * 

al/l  \m.  p 1 :•!/  [ VI . p\ 


If  the  next  approximationsdif f ered  from  the  previous  ones  by 
less  than  some  small  value  £.2,  the  solution  at  the  given  point 
was  taken  to  be  equal  to  the  latter  approximation.  Otherwise,  the 
iteration  process  was  continued.  It  is  necessary  to  note  that 
convergence  of  the  iteration  process  was  assured  everywhere,  except 


for  a section  of  the  time  axis  directly  adjoining  the  coordinate 


origin.  In  this  section  the  form  of  the  unknown  functions  was 

4 

given  beforehand  in  accordance  with  a well-known  approximate  solution 
As  an  example.  Fig.  20-  shows  the  results  of  solving  the  plane 
problem  for  (5^  = -0.1.  The  vector  D was  assigned  in  the  form 

0.0 
0,5 
1.0 
2,0 

4.0 

5.0 

12.0 

52,0  I 


Fig . 20  . 


r 


*** 

t 


f.  *. 
b ' ' . 


Consequently,  the  partition  points  concentrated  near  the 
coordinate  origin,  where  the  most  intensive  change  in  the  unknown 
functions  takes  place.  At  the  points  D[l]  ,D[2]  , D [ 3 ] the  values 
of  the  functions  were  assigned  from  the  approximate  solution, 
and  at  the  remaining  points  they  were  determined  using  the  described 
iteration  process. 

The  results  of  solving  individual  problems  using  this  program 
are  shown  as  graphs  together  with  the  results  of  calculations 
carried  out  using  approximation  methods.  Comparison  of  these 
results  is  made  in  the  paragraphs  devoted  to  the  corresponding 
approximation  method?. 


2.  Numerical  Solution  of  the  Problem  of  Heat  Exchange  During 
Flushing  of  a Well 


In  this  case  the  heat  problem  in  ground  also  reduces  to 
system  of  equations  (A. 11) - (A. 14) . Equation  (A. 15)  preserves 
its  form,  with  this  sole  difference:  in  place  of  T in  this  equation 

it  is  necessary  to  substitute  , the  temperature  of  the  flushing 
fluid  in  the  annular  space. 

The  basic  difference  will  consist  in  the  conditions  of  heat 
balance  in  the  well,  which  in  our  case  are  described  by  equations 
(3 . 89 ) - ( 3 . 90 ) . Let  us  transform  these  equations  in  the  following 
manner.  Differentiating  equation  (3.90)  with  respect  to  z,  we  have: 


4*Tt  (dr,  &Ti\ 

■3?)- 


(A. 17) 


125 


II 


whmbhhi i i 


Equation  (3.89),  taking  into  account  (3.90),  may  be  written  in 
the  form 


oT~. • <Uj 

■0T~.  a3l»+d:'  (A.18) 


Here 


2 r./.X, 
;G  • 


Substituting  (A. 18)  into  (A. 17),  we  obtain 

< r-T 


(A. 19) 


Thus,  adding  to  system  of  equations  (A . 11 ) - (A . 14 ) equations  (A. 18) 
and  (A. 19),  and  also  the  equation 


(A. 20) 


we  obtain  a complete  system  of  equations  for  determining  the  unknown 
functions.  This  system  was  solved  on  a BESM-3M  computer.  The 
program  was  written  in  ALGOL  for  the  c<  -translator . 

The  fundamental  ideas  on  which  the  program  worked  are  as  follows. 
The  entire  stratum  is  broken  up  into  q horizontal  sections.  On  the 
time  axis  n points  were  singled  out;  the  coordinate  of  the  mth  point 
was  denoted  by  D[mJ  , where  m = 1,  2,  ...,  n.  For  each  moment  in 
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time,  beginning  with  m = 2,  system  (A. 11) - (A. 14 ) , (A. 20)  was  solved 

for  all  the  horizontal  sections  using  the  method  presented  above. 
The  obtained  function  vQ  was  used  for  determining  the  functions 
T^  and  T ^ for  the  given  moment  in  time  from  equations  (A. 18) - (A. 19) 
The  function  T2  was  substituted  into  equation  (A.  20)-,  and  system 


(A . 11 ) - (A . 14 ) , (A. 20)  was  again  solved  for  the  next  moment  in  time. 


APPENDIX  2 


Program  fo.  Solving 

-(A.  16)  Using  the  Programming  Language  aluui; 

bf  a BESM-3M  Computer  ~ 


1.  begin  real  T.  a,  d,  f,  a I,  (?2,  f2,  f3,  e4, 

2.  nn,  vQ.  yl.  t'2,  b I,  b 2,  «2,  Rt,  Rnr, 

3.  integer  in,  n,  p,  q.  i,  j,  .V;  n?acLi".  q) : begin, 

4.  array  vOj,  ti /,  t'2/,  //,  s/,  7'g/,  lO/  1 , o I / 1 . :'2/l , 

5.  //I,  s/1  [ 1 : /!,  1 : g),  v.  1,  x2,  7.1,  7.2,  h,  FO f I : <7], 
ft.  el.  D[l  : «]; 

7.  real  procedure  /Ob/); 

S.  value  //;  real  ir. 

0.  begin  if  abs{ij)<.  10  then  go  to  \'l  else 

10.  .Ill:  begin  real  pc.  sum,  integer  r; 

11.  r:  = I ; su/;j;=0;  pc:  — — !/(SX'/l, 

12.  it:  sum:  = sm’f4-  ( — l)frX/'<-’;  r:  = f-f-I: 

13.  pc:  = pex  ( - (2X'  ~ I ) 1 2/ 1 SX'-X?/)  I ; 

14.  if  i7.bs  (pc)  > |i — 4 then  go  to  if  else 
15  /0;  ==  ( 1 -i-suin) Isqrt  (ahs(tj)  X6  2832) ; 

16.  go  to  P end  Aft; 

17.  A';:  begin  real  pc,  sum:  integer  r; 

18.  r:=l;  sum:  = 0;  /;<.•;=(;// 2)  f 2; 

19.  /./:  sum:  = sum-\-pc:  r:  = r-j-l; 

20.  pc:=pcy,  (ijl ( 2 X ' ) ) 1 2; 

21.  if  r<  10  then  go  to  l.l  else 

22.  10:  = (\—sum)Xcxp( — 11)  end  A'f;  P: 

23.  end. 

2 1.  real  procedure  / 1 (//) ; 

25.  value  //;  real  //; 

26.  begin  if  abs(i/)<10  then  go  to  .VI  else 

27.  .Ul:  begin  real  pc,  sum:  integer  r; 

23.  r:  = I ; su/ff.  =();  /Jtr:  = 3/ ( S X// ) : 

29.  iter:  sum:  = sum-ri  — 1)  f r'XPc;  r;  = r-f-  I; 

;0.  pc:  = pc  X (4  — (2X;‘ — 1 ) f 2)/(8XrX'f) : 

31.  if  ais ipe)  > 10"*  then  go  to  iter  else 

32.  /! ; = ( l+sum)lsqrt  (abs(tj)  X6.2832) ; 

33  go  to  PI  end  All; 


"'fM!  I-  ■ rn  t 

f viiitij  i t UlL  i 

• nnw  e * 
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31.  AM:  begin  real  pc,  sum;  integer  r; 
35.  /■:  = 0;  sum:  — 0;  pe:=y/2; 


33.  L 1 : sum:  = sum-\-lur,  r.  ==  r-f- ! ; 

37.  pc:=zpcX(yl 2)  f2/('X  (r+l)); 

38.  if  r<  10  then  go  to  LI  else 

39.  /I  :=5«mX«P( — !/)  end  'VI;  PI: 

40.  end; 

41.  real  procedure  /I):  value  x,  /l; 

42.  real  .r;  array  A;  begin  integer  j; 

43.  /:  = 1 ; TS:  if  *<D(/+l]  then  begin  W:  = 

44.  p]X.(x-D[j])+A[j,  p]X{D[j+\)—x))l 

45.  (D[/-H]— Df/]):  go  to  RS  end  else 

46.  /:  = / -ft;  go  to  TS;  RS:  end: 

47.  real  function  dlf.r,  //)  = (/Of  l/(2Xx 1 [P]  X(* — [)) ) ) — 

48.  /l(l/(2XxlMX(v— '/))))X'V/(i/.  vOi  1 ) ; 

•19.  real  function  d2(x,  u)  =exp\  — ( ’.V (;/.  sj)  — 1 ) f 2/(4 Xxl  [p]  X. 

so.  (jc— w)))x(/ o( vv’(;/.  s/)/(2Xxi[p]x<*— !/)))— «/)X 

51.  / ! ( IV' (//,  -•>7)/f2X>:l[P]X(.v-/y))nX^f.'/.  si)XW(y.  wljl): 

52.  real  function  dij  0(x.  y,  k)  = ( 1 / ( 2 X x.  I [p]  X (*—  y)  f2) ) X 

53.  (rfl(.r,  ij)—  2X(X— {/)X/l(l/(2Xxl[p]X(-v— !/>>)X 

X ( vv ' ( y 4-  e . 

54.  fj)-Vi'(y,fj))/K-d2(x.  *t)); 

55.  real  function  d3(.v,  u)=exp( — (tt  (.v,  s/)  — 1 ) f 2/ ( 4 X 

56.  xl[p]X(*— 0)>)X(^(*.  sj)Xl 0(U"(Jf.  s/)/(2Xxl[/>]X 
X (x—t/))) 

57.  --/ 1 < sj)H2X*  1 [p]  X (Jf— £/) ) ) f X ({/.  eO/1)  ; 

58.  real  function  i/4(.v,  t/)  =exp(  — ( Uvf.v.  Sj)  — W(y , s/))f2 

59.  /(4Xx  I [p]  X (-^ — i/) ) ) / (2Xx  l [p]  X (*—</)  1 2)  X 

(. — / o ( ir  (x,  si)x 

60.  W( y,  s / ) / ( 2 X x 1 [p]X(x~y)))XW(x,  sj)+mW(x,  sj)X 

61.  W, y,  s/)/(2X>'-l [/’] X (x — y)))XW (y,  s/))X »'(!/.  «/)X 

62.  U"f(/.  wl/'l); 

63.  real  junction  d/7l(.r,  y,  r)  = (l/(2Xxl  [p]  X (* — y)  f2) ) 

64.  Xd3(x,  y)—exp(—(\V(x,  sj)  — 1)  f2/(4Xxl  [p]  X 

55.  ( x-y) ) ) / ( y.  1 [/;]  X (x-y) ) X/ 1 ( W(x.  sj)f(2X*  1 [p]  X 
X ( x—y) ) ) 

66.  X(H7(|/+r..  fi)-W(y,  fj))h+dl(x,  y); 

67.  real  function  dif 2(x.  u)  = \V(y,  sj)X'X,(i/,  r.’ 2 / 1 ) X 

68.  exp  (-  ( IP  (x,  sj ) - \lv  (y.  sj) ) 12/(4 X x2  [p]  X ( x-y ) ) ) X 

69.  ( U'<r,  sj)XfO(W(x,  sj)XW(‘J,  w)/(2Xv.2f/;J  X (•«—?/)))  -- 

70.  nv(H.  s/)xn(W'(.v,  *j)XW(y,  s/)/(2Xx2[p]X(*-v))))/ 

71.  (2Xx2[p]X  (•'•-'/)  t2); 

72.  real  function  z0(yj  = dif0(D[/u  , u,  I0"2); 

73.  real  function  *1  (y)  =dif\  (D  m , ?/.  1 0 2) ; 

74.  real  function  z2(y)  =dif2(D  in  , </); 
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" rial  function  /'(.v)  =-  (7)1  X{*X(1-  <x2,‘2) -r«2X*X 
7u  !!i(.c))-|-W/.r)/(x2[/»jX^im]J; 

77.  real  function  c.’J  ( f/ ) ~ijY.l (y)  X?XP{-  W(m,  nj  — 

78.  !ni->H-\X*2\p]XD[ir.]))+n  (sj[m,  p]XyH  2\ 

79.  v:2[p\XD[m\))\ 

80.  reaJ\T.  n,  xl.  x2.  /.I.  72,  d.  F0,  f,al,o2.e\,  f2, 

81.  //,  r3,  D,  ].  1 1,  //.*;,  Rm,  a2,  Tgj,  si,  vOil,  t’l/I,  i»2/l); 

82.  for  ,o.  = l slep  ! until  q do 

S3.  BO:  begin  Rt:  = 4Xx2[p] X (a2Xin (Rw)  -r 1 ) \2/(RmX 
84.  (2X (a2Xln {Rm) -l){ 2-2Xa2X («2Xln t Rm)  + 1 ) -r 
-~2{2X 

83.  (2 — 2 X v.2 — -st2 1 2 > / (Rm)  f2)4): 

86.  b I : = a2+«rx  ( Tgj  [I,  p\ -FG[p]  )U27<y.2\p]  X«»«X(*2X 

87.  In  (Rm)  -f  M , 2) ; 

88.  b2.='x2-\-R:y,{Tgj[\.  p]  — F0[p])X ( (2X^2— a2f2— 2)/ 

89.  Pm  -cc2  X&'n  X ft  - a2+ a2X  In  t Rm) ) W ( 4 Xx2  [ p]  X 

99.  (s2Xlni/?»i)  + l)f3)+*2X(/r0[p]~r^[l1  ;;])/(a2X 

91.  In  (Rm)  -4-1); 

92.  tO/(I,/»  = vO/1  f I , pi;  cl/[  l,  />]:  = cl  ;i  [ 1,  p]; 

93.  ; 2/[  I . />j:  = c2/i  [I,  nj  end  Z?0; 

94.  Rm;  begin  for/;=l  step  1 until  q do 

95.  R 2:  begin  for  m:  — j step  ! until  n do 

96.  7?4 ; begin  Z'2j[j — 1,  />]:  = c2/l  [/ — I.  p]\  t’i/[/—  I.  p] 

97.  : = i’lil  [ j — 1,  /;) : 

98.  i<0'.  = sini/>s(j0,  t),  Dfui  — 1],  t4,  1 4) ; 

92.  v \:--=  simps  (z\,  0.  D|/n— 1 J . el.  f.4); 

100.  u2  = si:/:/) «(<•_>.  0.  /;[/»  — !],  f4,  f4); 

!0I . U:  ;jJ:  = i-0+sj/n/)Sfr0,  D\m—  1].  — /. 

1 02.  f.5.  f3)  + sirups  U0,  D(n]—f,  D( vi J — r 1 [■•«] , f3,  f3); 

103.  c0/ [w/.  r»] - — 1-0/ [/n,  p]-4-H  f'm]/(s<?rf(2»—  l)X<(s‘/rt 

101.  1 2)— 2)'XMi /.>[;/(]  — f1  [m])-f;0(D ["i]—; k1  [m},2> ) — 

1 0.5.  fj[l , p ] / ( x 1 [p | XD [m]  | X/U  1 !(y- 1 [/’]  XV \ '» j X2) ) ; 

106.  z‘lj[/n,  p]:=zvl+siinps(zl,  D\m—  ij,  D[:n]—f, 

107.  e3,  , 3)  -r-simps  XI,  77) [m]  — f.  D [o;] — f 1 [:«J , 

108.  t'3.  f.3);  l /[«!,  p] -f  f.I  f"']/t 

1 09.  sqrt  ( 2)  - - 1 ) X ( l (2 ) -2)  X*  1 iD  l <»  J ~ f 1 [«« ] ) 4- 

i 110.  zl[/«]— f I fm]/2) ) — fj[l.  p]Xcxp(—(si[m,  p}  — 1 ) f 

111.  2/(4X-/.lfp]XO[«i]))X/I(i/K  p]/(2X/-l[pjX 

112.  DM))/ixif,v]x£>[''i]); 

: 113.  fj\\m,  /)]:=(  1/a) Xi'0/[//».  pH-T'g/r.n.  p], 

114.  *c2;|»»,  p].  =i'2-fsimp.s(z2,  D[/:t— !].  D[m|— /, 

115.  f3.  f.3)  -{-simps  (c2,  D[w]— /,  D\m]~  f!  [ml. 

116.  . 3,  k3)+f1  [n: ]/(.'C,’r.' (2)  — 1 ) X f {sqrt  {2) — 2)X 

117.  z2 (!) { m] — k 1 [ m j)  -v-<2(D[m j — r l [m]  12) ) + simps 
1 IS  i <3.  1.0,  Rm,  f3,  f3); 
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119.  s/I  [//i,  p\:  =ij[m- ••  I,  p]  \ (I) [m  | — D[; n-  - 1 ] )/2X 

120.  (— X/[/JjX(i’l/f"1.  /’]  + - I 1,  p])-H2[p]X 

121.  (1-2/ 1 mi,  p +i'2/[/n-  -1.  />])); 

122.  vO/l  [mi,  /)J.  =0.5X,(^'0/1  [mi,  pJ+^'O/Tmi,  pi); 

123.  vl/1  ,vi,  p :=s0.5X(l’l/l  mi,  p -f t'l/  mi,  p ); 

124.  i*2/l  mi,  p]  =0.5X  (u2/l  [mi,  p]  -\-v2j  m,  p ]); 

123.  //I mi,  pJ:  = 0.5X <//['«,  PJ+//I  [mi,  p]); 

120.  a'/ [mi,  pJ;  = 0.5X(^/[Mi,  pj+s/’l  [mi,  p] ) ; 

127.  torite.  (i'0/I  [mi,  p],  i'1/1  [m,  p],  v‘2jl  [n,  p], 

128.  Tgj[m,  p],  ij[m,  p],  sj[m,  p]); 

129.  D:  il  al>s(v\  j.\[m,  p] — v\ /[mi,  p] ) -\-abs  [vOi\ 

130.  [mi,  p] — i'0/[mi,  p])>e2  then  go  to  ds 

131.  else  end; 

132.  for  mi:=!  step  1 until  n do  begin 

133.  Tgj[m,  p):  = Tgj[m , p]  + (al  Xw0/'[mi,  p]—  fl2) 

131.  end; 

135.  1 0/ [ I . p ]•'  = m 1 / [ 1 . p]:=—c/.XTgj[\,  pj; 

136.  i’0/1  [1,  p ] : = v I / 1 [I,  p]:==;>0/[l,  p] ; end  R 2 

137.  end  go  to  Rot  end  end 
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